One Step at a Time: How Model Time Step Significantly Affects Convection‐Permitting Simulations

Abstract

We show that there is a strong sensitivity of cloud microphysics to model time step in idealized convection‐permitting simulations using the COnsortium for Small‐scale MOdeling model. Specifically, we found a 53% reduction in precipitation when the time step is increased from 1 to 15 s, changes to the location of precipitation and hail reaching the surface, and changes to the vertical distribution of hydrometeors. The effect of cloud condensation nuclei perturbations on precipitation also changes both magnitude and sign with the changing model time step. The sensitivity arises because of the numerical implementation of processes in the model, specifically the so‐called “splitting” of the dynamics (e.g., advection and diffusion) and the parameterized physics (e.g., microphysics scheme). Calculating one step at a time (sequential‐update splitting) gives a significant time step dependence because large supersaturation with respect to liquid is generated in updraft regions, which strongly affect parameterized microphysical process rates—in particular, ice nucleation. In comparison, calculating both dynamics and microphysics using the same inputs of temperature and water vapor (hybrid parallel splitting) or adding an additional saturation adjustment within the dynamics reduces the time step sensitivity of surface precipitation by limiting the supersaturation seen by the microphysics, although sensitivity to time step remains for some processes.

Introduction

Convection‐permitting models are used in a large number of different meteorological applications and have provided a step change in rainfall forecasting (Clark et al., ) over previously used models with parameterized convection. The convection‐permitting models are commonly used for numerical weather prediction (e.g., Baldauf et al., ; Tang et al., ), examining severe weather case studies (e.g., Hanley et al., ; Kunz et al., ; Warren et al., ), process studies (e.g., Barrett et al., ; Barthlott et al., ), and to assess the changes of precipitation under certain different external forcing, such as aerosol influences (e.g., Fan et al., ; Rieger et al., ) and climate change (e.g., Kendon et al., ). Because of the widespread use of convection‐permitting models it is therefore very important that the relevant processes within the model are accurately represented by the explicit and parameterized physical equations, which includes the correct numerical implementation.

The current suite of convection‐permitting models are more skilful at predicting precipitation than coarser resolution models with parameterized convection (Clark et al., ; Kendon et al., ; Mass et al., ; Mittermaier et al., ; Roberts et al., ; Roberts & Lean, ). However, there are still a number of weaknesses that are evident in the models, particularly in their representation of convective precipitation structures and the correct timing of the onset of convection.

There are numerous potential causes for these weaknesses, such as convective initiation affected by interactions with the land surface (Hanley et al., ; Schneider et al., ) or errors in simulation at the synoptic scale (Barrett et al., ), only partially resolving the convective structures because of model resolution (e.g., Chemel et al., ), the weakness of single‐moment microphysical parameterizations compared to double‐moment schemes (Igel et al., ), or the inherent unpredictability of convective events (Barthlott et al., ).

Previous studies have addressed how the model skill and representation of cloud systems change as the horizontal resolution (Chemel et al., ; Hanley et al., ; Stein et al., ) or vertical resolution (Barrett et al., ) are changed. However, to our knowledge, no study has so far looked at the impact of changing the model time step independently from resolution changes. In this paper, we document the significant sensitivity of precipitation to the chosen model integration time step and in which way the model numerics, in particular the operator splitting, contribute to this sensitivity.

The operator splitting is commonly used to split fast modes (e.g., gravity waves and sound waves) from the slow modes (e.g., advection) within the dynamical core of the model (Burridge, ; Gadd, ; LeVeque & Oliger, ). The fast and slow modes are each solved separately. This is advantageous because the fast modes require a small time step to ensure accuracy and stability, whereas the slow modes can be accurately simulated with longer time steps. The splitting of the fast and slow modes is therefore computationally efficient. The most common way to implement such splitting in the model is to iterate the solution for one of the two modes first. After this first mode has been solved for, the second mode is calculated using the updated fields as input for the calculation. Strictly speaking this method is only valid if the two subprocesses are commutable, meaning that the solution is identical regardless of the order of the two processes. However, accurate solutions of the Navier‐Stokes equations are nevertheless achieved using such splitting of fast and slow modes (e.g., Skamarock & Klemp, ).

Within weather prediction models, the resolved fluid mechanics representing the advection and diffusion processes (“dynamics”) are solved separately from the parameterized physical processes such as cloud microphysics, radiative transfer, and turbulent boundary layer processes (“physics”). The numerical method of solving these two separate parts often involves “operator splitting” (sometimes referred to as “additive splitting,” e.g., in the Weather Research and Forecasting, WRF, model documentation; Skamarock & Klemp, ; Skamarock et al., ); however, we use the term sequential‐update splitting following the terminology of Donahue and Caldwell (). In many cases, the parameterized physics provide tendencies that are given to the dynamics part of the model as forcing terms. However, the cloud microphysics is sometimes treated differently from other parameterized physics, where it is calculated after the completion of the rest of the dynamics and physics calculations (e.g., in the WRF model, Skamarock & Klemp, , and in the COnsortium for Small‐scale MOdeling, COSMO, model, Baldauf et al., ). This allows the microphysics to react to supersaturation produced or cloud water condensed during the time step. However, it also means that the input values for the microphysics parameterization are dependent on the output of the dynamics and physics parameterizations. In convection‐permitting simulations, which resolve large vertical velocities and mean that significant supersaturation with respect to liquid (hereafter SSLiq) can be generated through adiabatic cooling in updraft regions within a time step, these supersaturations are then passed to the microphysics scheme. A flow chart of the sequential‐update splitting used in the COSMO model is shown in the left column of Figure .

Flow chart illustrating the numerical implementations of the three different splitting types. The flow of the model variables (in blue) temperature (T), water vapor mixing ratio (qv), and a dummy variable representing all hydrometeor quantities (qx) through the different processes in the model is shown. Processes written in capitals change the value of the model variables. The latent heating illustrated only reflects the heating due to phase changes in the microphysics parameterization. Additional latent heating from condensation is treated within the saturation adjustment.

In some models, SSLiq generated in updrafts is not passed to the microphysics; instead, all SSLiq is condensed to liquid water via a “saturation adjustment” (Lebo et al., ) before the microphysical processes are calculated. In these cases, SSLiq in the cloudy regions is relatively insensitive to changes in the model time step; however, supersaturation with respect to ice (hereafter SSIce) does still exist and does still depend on the model time step in regions where SSIce > 0 but SSLiq < 0. The amount of condensation produced via saturation adjustment also scales linearly with the time step length. A modified version of COSMO, where saturation adjustment occurs within the model dynamics (and therefore before the microphysics call), is illustrated in the middle column of Figure . Some other models use a subgrid cloud scheme, where part of the grid box can contain cloud while the other part remains cloud free—in these models the condensation parameterizations typically assume instantaneous removal of SSLiq and, in this sense, are similar to saturation adjustment.

An alternative to the sequential‐update splitting described in the above two paragraphs is parallel splitting. With parallel splitting, the fast and slow modes are treated separately and independently, and the new solution is found by linearly adding the contribution of the two modes together. Therefore, there is no specified order to these processes, but rather, they are deemed to occur at the same time. Parallel splitting is only rarely used because sequential‐update splitting is often found to be more numerically stable and give better results. For example, Rood () found that parallel splitting introduces small numerical errors in chemistry transport models but that the errors are much smaller than errors coming from inexact advection. Williamson () and Beljaars et al. () also compare parallel and sequential‐update splitting techniques and find that sequential‐update splitting gives a better solution. Beljaars et al. () further advocate the use of sequential‐tendency splitting. Dubal et al. () show that the parallel splitting maintains a steady state solution regardless of time step when explicit time stepping is used, whereas the sequential‐update splitting does not; however, parallel splitting gives a biased solution with implicit time stepping. More recently, Donahue and Caldwell () looked at the importance of the order of processes in a climate model. Although they found that the model is rather sensitive to the order of processes when sequential‐update splitting was used, parallel splitting was not recommended either due to problems “keeping concentrations positive.” However, none of these studies investigated models where convection was explicitly resolved and therefore the contributions from vertical motion and the resulting condensation have been overlooked to date.

This paper is organized with the description of the numerical model setup in section . The sensitivity to model time step is shown and explained in section for COSMO simulations, and for a simple differential equation in section . Section also shows that using parallel splitting gives better results when solving the differential equation numerically. This method is then applied to COSMO simulations, with the results shown in section . Analysis of simulations with an additional saturation adjustment during the model dynamics is also shown in section . Section looks at the causes of the sensitivity and assesses the changes to microphysical process rates. Section further discusses the significance of the time step sensitivity. Section contains conclusions.

Numerical Model Setup

The COSMO model (Baldauf et al., ) is a nonhydrostatic model, developed and used operationally by the German Weather Service and other European national weather services for forecasts on limited area domains. For simulations in this paper, COSMO version 5.3 has been used in an idealized setup, with 1‐km horizontal grid spacing. There are 64 vertical levels giving a vertical grid spacing of 300 m at 5‐km altitude and 400 m at 10‐km altitude. There are open boundary conditions, no surface fluxes, and no radiation. The initial wind and temperature profiles are taken from Weisman and Klemp () with a near‐surface moisture content of 12 g/kg and unidirectional wind shear. The winds are westerly and increase from 0 at the surface to 31 m/s at and above 6‐km altitude. The convection is initiated with the release of a warm bubble at the beginning of the simulations, with horizontal radius 10 km, vertical radius 1.4 km, and a temperature perturbation of 2 K in the center of the bubble decreasing to zero at the edge, following a cos‐squared distribution in space.

Operator Splitting

The default order of the model calculations follows sequential‐update splitting (as illustrated in Figure , left column); therefore, the model dynamics (advection and diffusion) is calculated first, the model variables are then updated, and these newly updated model variables are used as input for the microphysical parameterizations. At the end of the microphysics, any remaining SSLiq is condensed to cloud water through the saturation adjustment. When using the two‐moment scheme in COSMO, the default setup is to use the saturation adjustment only after the microphysics.

Two modified versions of the model numerics are tested in this paper. The first of these is sequential‐update splitting with an extra saturation adjustment during the model dynamics as part of the Runge‐Kutta stepping (hereafter SUS+satad; Figure , middle column). Simulations where the saturation adjustment is included in the model dynamics are discussed in section . The second is a hybrid implementation of parallel splitting, which only considers the microphysical parameterization and the model dynamics. Our implementation of parallel splitting in COSMO (right column of Figure ) is not strictly parallel but aims to avoid potential numerical issues where the cloud variables could become negative if they are simultaneously advected out of a grid box by the dynamics and also depleted by the microphysics parameterization. Therefore, the dynamics is indeed calculated first, but the initial values of temperature and water vapor mixing ratio for the dynamics are also used as the initial values for the microphysics. The hydrometeor mass and number densities are advected by the dynamics and then passed directly into the microphysics. The tendencies of temperature and water vapor mixing ratio from the dynamics (“dynamics increments”) are saved, but not applied to the prognostic variables until after the microphysics parameterization is completed. Simulations using this parallel splitting are discussed in section .

Cloud Microphysics

Simulations shown in this paper use either the single‐moment microphysics scheme from the operational model (Baldauf et al., , section 2b) or the two‐moment, six‐category (cloud, rain, ice, snow, graupel, and hail) microphysics scheme (Seifert & Beheng, ) including the addition of the hail category from Blahak (). The two‐moment scheme predicts number and mass concentrations for each of the six microphysical categories, and each use a generalized gamma distribution to describe the shape of the particle size distribution. The microphysics parameterization has, by default, the same time step as the native model time step. A short description of the processes relevant to discussion of the time step dependence are included below, ordered by the order in which they are calculated within a time step. For a more complete overview of the microphysics scheme, readers are referred to Seifert and Beheng ().

Order of Processes

Each of the microphysical processes is calculated using sequential‐update splitting. The calculation of SSLiq and SSIce occurs at the beginning of the microphysics parameterization. The ice phase processes are calculated first (nucleation, freezing and vapor deposition, and collection) followed by the mixed‐phase processes (e.g., riming and melting) then melting and evaporation. After that, the liquid phase processes (autoconversion, accretion, collection, breakup, and evaporation) are calculated followed by saturation adjustment and sedimentation. A flowchart summarizing the order of processes can be seen in Figure 1 of Seifert and Beheng ().

Condensation via Saturation Adjustment

Condensation is handled within the saturation adjustment part of the model; there is no explicit calculation of the condensation rate. In the default COSMO setup using the two‐moment scheme, the saturation adjustment is only active after the microphysics calculations have been completed. In the SUS+satad setup, an additional saturation adjustment takes place during the model dynamics. The latent heat changes related to condensation are also applied by the saturation adjustment and not where the rest of the latent heating is applied after the microphysics.

CCN Activation

The cloud condensation nuclei (CCN) activation uses the Segal and Khain () lookup tables, which determines the number of activated cloud droplets based on the vertical velocity. Activation of CCN occurs only where SSLiq ≥ 0 and grid‐scale vertical velocity is positive, although it is the vertical velocity, rather than supersaturation, that is used to determine the number of CCN that are activated. A term restricting CCN activation to grid boxes where SSLiq was increasing with height was removed because it often resulted in activation occurring in alternate grid boxes with height. Activation of CCN can occur both at cloud base, and within the cloud at grid boxes where SSLiq exists at the beginning of the cloud microphysics calculations (see Figure ). Default CCN concentrations in this paper are “maritime” conditions, with 100 CCN cm⁻³. The CCN concentrations are constant in the lowest 2 km, decreasing exponentially with height above 2 km, with a scale height of 2 km. The sensitivity of increasing the maximum CCN concentration to 1,700 cm⁻³ is also shown.

Ice Nucleation

The ice particle number concentration (n_ice; equation ) is parameterized following Fletcher () with a modification from Huffman and Vali () to account for the effect of vapor depletion in the vicinity in the presence of other ice particles. [Image Omitted. See PDF]

In this equation, N₀ = 10⁻² m⁻³, S_i is the fractional supersaturation with respect to ice, S₀ is the fraction supersaturation with respect to ice at water saturation, B = 4.5, k_f = 0.6 K⁻¹, and T is the air temperature in kelvin. If B = 0, then is the Fletcher () parameterization. Ice multiplication through the Hallett and Mossop () process is also parameterized.

Riming

The collection of liquid particles by falling ice particles is parameterized by a sweep‐out kernel and a collection efficiency. The sweep‐out kernel calculates how many collisions between an ice hydrometeor species and liquid hydrometeor species are expected using the horizontal cross section of the particles and the difference in their terminal fall velocity. This is integrated across the particle size distributions of both the ice and liquid hydrometeor species in question. Not all of these collisions will result in riming, the fraction that do is parameterized by the collection efficiency E_coll: [Image Omitted. See PDF] where ${\overline{D}}_c$ is the mean liquid‐particle diameter, ${\overline{D}}_0=15$ μm, and ${\overline{D}}_1=40$ μm. The scheme neglects differences in the collection efficiencies of number and mass. The complete equations are given in Seifert & Beheng (, equations 61–67).

Autoconversion and Accretion

The autoconversion of cloud drops to rain and the accretion of cloud drops by falling rains drops are parameterized following Seifert and Beheng (). The autoconversion rate approximately scales with $L_c^2{\overline{x}}_c^2$, where L_c is the cloud water content in kilograms per cubic meter and x_c is the mean mass of cloud drops. The full equations are given by Seifert & Beheng (, equations 4–6). Accretion rate approximately scales with L_cL_r, where L_r is the rain water content in kilograms per cubic meter. The full equations are given by Seifert & Beheng (, equations 7 and 8).

Clipping of Values

Within each of the microphysical processes, the tendencies are limited to the maximum available mass and number of the relevant hydrometeor category. If the parameterized conversion rate exceeds the available mass or number, the rate is reduced. The mass and number densities are forced to be nonnegative both before the microphysical calculations and after the microphysical calculations and sedimentation are complete. Furthermore, minimum and maximum values for mean particle mass are specified. The number density is adjusted if the mean particle mass falls outside this range.

Time Step Dependence in COSMO Simulations

The effect of the model time step is demonstrated by using the idealized model setup and varying the model time step from 1 to 15 s, whereas all other model parameters are fixed. This range is chosen because excessive numerical diffusion becomes an issue for time steps shorter than 1 s at this horizontal resolution, and numerical stability becomes an issue for longer time steps. The total precipitation (rain, snow, graupel, and hail, although no snow or graupel reaches the surface) and hail reaching the ground in each simulation are shown in Figure . These plots reveal three undesirable characteristics:

1. The total amount of precipitation and hail decreases significantly as the time step becomes longer.
2. The location of the heaviest precipitation and greatest hail accumulation varies strongly and systematically as the model time step is changed.
3. The effect of increasing CCN number concentration varies depending on the chosen time step. Not only the magnitude but also the sign of the aerosol‐cloud interaction effect changes.

Summary of precipitation and hail after 2 hr of simulation time, for simulations with different time steps. Upper panels: Two‐dimensional x‐y maps of total precipitation (top row) and hail fall (second row) for simulations with 100 CCN cm−3. The domain shown covers 140 km in the east‐west direction and 36 km in the north‐south direction. Lower panel: Total precipitation and hail summed over the model domain. The arrows show the effect of increasing CCN concentrations to 1,700 cm−3. CCN = cloud condensation nuclei.

The surface precipitation total decreases by 53% (Figure , top row) and surface hail almost vanishes (Figure , second row) as the time step length is increased from 1 to 15 s. Additionally, the maximum values at any individual grid point are also strongly affected by the model time step. The maximum precipitation decreases from (25.1/20.8 mm) to (12.0/8.0 mm) for (clean/polluted) simulations when the time step is changed from 1 to 15 s.

The location of the precipitation maximum is also strongly affected. In the 1‐s time step simulation, the largest precipitation total are directly underneath the two convective updrafts; this is also the region with the most hail. In the 15‐s time step simulation, the precipitation maximum remains along the central axis of the model domain, and there is little‐to‐no hail reaching the surface.

For a model time step of 1 and 2 s, there is more precipitation in the cleaner (100 CCN cm⁻³) environment. With longer than a 2‐s time step, there is almost no sensitivity of total precipitation to changing CCN. For time steps between 3 and 6 s, there is slightly more precipitation in the more polluted environment (1,700 CCN cm⁻³). A similar swap of the sign of the aerosol impact is seen for hail. For time steps shorter than 3 s there is more hail in the cleaner environment, for medium time steps there is more hail in the polluted environment, and for 10 s and longer there is almost no hail in either simulation.

The effect of the time step on these three aspects of the model output may significantly contribute to the continuing challenges of (1) quantitative precipitation forecasting, (2) pinpointing the locations of extreme precipitation, and (3) deducing the aerosol‐cloud interactions in convection. Furthermore, the microphysical behavior of the simulated cloud also varies as the time step is changed.

In addition to the surface precipitation, the vertical distribution of different hydrometeor types can be largely affected by different time steps. Figure shows an example of the vertical distribution of each of the six hydrometeor classes in COSMO for short (1 s) and long (15 s) time step simulations. In the short time step simulations, there is significant supercooled liquid cloud water and rain water. However, in the long time step simulations, much of this supercooled liquid in the cloud is frozen. The difference is particularly clear between 4‐ and 7‐km altitude. As a result of this freezing and the associated additional latent heat release, the storm updraft was strengthened (not shown) and the simulated cloud top height increased (Figure ).

The vertical distribution of hydrometeors after 60 min of simulation time. (a) Simulations with a time step of 1 s; (b) 15‐s time step. The colors show the contributions of the each hydrometeor category.

This substantial change to what is happening within the cloud reveals that there is a large, systematic problem with the long time step simulation in the default setup of this model that affects both typical model output variables (such as precipitation and hail) and also significantly affects microphysical processes and, thus, how the simulated storm evolves.

Substepping the microphysics, by calculating the process rates and sedimentation using numerous shorter time steps within a single dynamical time step, does not improve the time step sensitivity. This is because the time step sensitivity is primarily caused by the model numerics and the input values for the microphysical parameterization, and not by the microphysical processes themselves. This will be explained in more detail in sections , and .

Time Step Dependence in a Simplified Model

The cause of the sensitivity to the model time step relates to the numerical implementation of the model, in particular to the sequential‐update splitting applied within the model. Sequential‐update splitting means that within one time step the model “dynamics” (the iteration of the Navies‐Stokes equations; mainly incorporating the advection and diffusion within one time step) is calculated first. After the dynamics, the prognostic variables are updated. The newly updated prognostic variables are then used as the input to the microphysics parameterization.

The reason that this leads to a large time step sensitivity in COSMO is explained in section . First, we will demonstrate with a very much simpler model for cloud microphysics why this causes a dependence on the time step. The simple model is shown in . [Image Omitted. See PDF]

The single equation represents the change in cloud water through two processes. First, the increase in cloud water with time from condensation (C; a fixed condensation rate of 0.002 kg·kg⁻¹ ·s⁻¹) caused by, for example, adiabatic cooling in a cloud updraft. Second, the conversion of cloud water to rain water through autoconversion—which is nonlinearly dependent on the current cloud water content—as often parameterized in cloud microphysics schemes (e.g., Khairoutdinov & Kogan, ; Seifert & Beheng, ). We use a value of A = 0.5. Once the cloud water is converted to rain water, it is removed from this simple system.

We implemented the sequential‐update splitting numerical method for solving —where the condensation term is calculated first, the cloud water content updated, and then the autoconversion term is calculated using this updated value. The model is then solved using several different time steps that range from 1 to 10 s. The output from these simulations (Figure a) reveal that there is a significant time step dependence in this simple system.

Time series of the evolution of cloud water content for different time steps in two different model setups. (a) Sequential‐update splitting. (b) Parallel splitting. The model time step ranges between 1 and 10 s, as shown in the legend. The analytical solution to is shown in both plots by the black dashed line.

As the model time step is increased, it is clear from the time series plots (Figure a) that the equilibrium value is different based on the time step used. For longer time steps the equilibrium value is systematically decreased.

The reason for this can clearly be seen when the numerical scheme is analyzed. The equilibrium value for the analytical solution to is $\sqrt{C/A}$. However, the equilibrium value for the sequential‐update splitting numerical scheme is $\sqrt{C/A}-C\mathrm{\Delta }t$. This difference occurs because the amount of condensation calculated in the first step of the sequential‐update splitting increases linearly with the longer time step. As a result, the intermediate values of q_c depend on the time step. Hence, the calculated autoconversion rate is larger when the time step is longer because the q_c input is larger. The time step sensitivity is then exacerbated by applying this increased autoconversion rate at a constant rate over the whole long time step.

An alternative methodology to solving such a system numerically is parallel splitting. This is similar to sequential‐update splitting, except that the input values for each of the processes is the same, and the tendencies from each process are then added together when updating the values of the prognostic variables at the end of the time step. This means calculating the tendency of cloud water due to condensation and the tendency of cloud water due to autoconversion, and adding the results to update the cloud water value. Figure b shows that when using this parallel splitting for the simple model that the numerical solution is consistent for all model time steps chosen. Furthermore, the value at equilibrium is identical to that of the analytical solution for ( $=\sqrt{C/A}$). The solution is the same for all model time steps because the input to each part of the model is consistent and does not depend on the time step. The parallel splitting method is able to represent the balance of these two processes in the simplified microphysics model. In the next section we show the impact of using this splitting method in the COSMO model.

COSMO Simulations With Modified Numerics

In this section, simulations with different numerics are evaluated. Two numerical alterations are made—the first is to include saturation adjustment within the model dynamics (middle column in Figure ). The second changes the setup to parallel splitting (right column in Figure ).

Additional Saturation Adjustment Within Model Dynamics (SAS+satad)

An alternative to passing large SSLiq values directly to the microphysics, is to perform an additional saturation adjustment within the model dynamics (the results presented above include saturation adjustment only after the microphysics). With the additional saturation adjustment, all SSLiq is removed, leaving previously liquid‐supersaturated air exactly at liquid saturation. The excess water vapor is converted directly to cloud water, and the temperature is adjusted to account for latent heating. This is usually argued to be a good approximation, because the timescale to condense liquid water is short in comparison to most model time steps. With saturation adjustment, the problem of passing large SSLiq values, with the amount being proportional to the model time step, to the microphysical scheme is removed and replaced by a similar problem: the passing large cloud water content values, which are also proportional to the model time step. This setup introduces other time step‐dependent effects and is numerically very similar to the simple condensation‐autoconversion model introduced in section .

However, COSMO simulations with the additional saturation adjustment within the model dynamics are less sensitive to the model time step than simulations with the saturation adjustment only after the microphysics (bars in lower panel of Figure ). However, the precipitation location and intensity (top row of Figure ) and the sign of the sensitivity to aerosol variations (arrows in lower panel of Figure ) remain sensitive to the chosen model time step.

Same as Figure , but for the sequential‐update splitting method including saturation adjustment within the model dynamics.

Parallel Splitting

Following the discussion in section , the sequential‐update splitting in COSMO was replaced by parallel splitting, but only for water vapor and temperature in the coupling between microphysics and dynamics. This means that the input temperature and water vapor concentration to the microphysics scheme is the same as the input to the model dynamics calculations.

This was achieved technically by saving the temperature and water vapor mixing ratio variables used as input for the dynamics and also saving the tendencies produced by the dynamics. The microphysics parameterization is called with the same input values as the dynamics and tendencies from the dynamics are added immediately after the microphysics calculations are finished (illustrated in Figure , right column). The main effect of the modification is that the calculated supersaturation depends much less on the model time step. The hydrometeor mass and number mixing ratios are advected by the dynamics before the microphysical calculations, which avoids potential problems of numerical instability from negative mixing ratios.

Switching to this parallel splitting implementation gives a much reduced time step sensitivity (Figure ). The three undesirable features of the default model setup have now been largely or completely removed for total surface precipitation, which varies by only 20% for all time steps used. Also, the systematic increase toward the shortest time steps has been reduced. The location of the heaviest precipitation is always within the two, splitting convective cells, which means that the precipitation distributions are now much more similar for all time steps. The impact of changed aerosol concentrations has a consistent sign and magnitude with the cleanest conditions resulting in 5–20% more precipitation with the precipitation maximum increasing by 10–40%. For hail, the total is 50–108% higher in the cleanest conditions, with the maximum at a single grid point higher by between 17% and 70%.

Same as Figure , but for the parallel splitting method.

Therefore, the use of parallel splitting is a significant improvement for consistency of results across all model time steps. However, there still remains some sensitivity to the amount of hail that reaches the ground; significantly more hail occurs for the shortest model time steps. This increase is also responsible for the increase of total precipitation at the shortest time steps (Figure ).

The reason for the remaining time step sensitivity of hail reaching the ground is likely that the individual processes in the microphysics parameterization are still implemented using sequential‐update splitting (see Figure 1 of Seifert & Beheng, ). This means that each process takes the mass contents of the hydrometeor species, calculates the changes caused by the currently considered process, and updates the mass (and number) contents accordingly. The next process then acts upon these updated values and further changes these values based on the new process. All of the numerous processes (e.g., 14 separate processes, each involving multiple hydrometeor species in Seifert & Beheng, ) can and will therefore interact with each other differently depending on the time step. An additional consequence of this is that the order of the microphysical processes within the microphysics scheme is important. Changing the order of the processes will change the input values for the processes, which can significantly change the resulting calculations as already demonstrated above. For example, if a process is very dependent on the presence of liquid water (e.g., riming or autoconversion), whether this process is calculated before or after the process of cloud droplet freezing and/or ice nucleation is decisive.

Causes of Time Step Dependence

In this section we discuss the causes of the time step dependence. First, by explaining how the time step sensitivity arises and, second, demonstrating how the microphysical process rates are affected.

Physical Interpretation of the Time Step Dependence

The time step dependence occurs because of how the model solves the physical equations in convective updraft environments. In the updraft regions, the adiabatic cooling is significant. A 10‐m/s updraft, for example, relates to a dry adiabatic cooling rate of 0.1 K/s, and model updrafts can be as strong as 50 m/s. The rapidly rising, and therefore cooling, air within the updraft leads to the production of large SSLiq. If the prognostic variables (most importantly temperature and water vapor mixing ratio) are updated first before the microphysics parameterization is called (as with sequential‐update splitting), then these large SSLiq values are passed directly to the microphysics parameterization. The SSLiq values passed to the microphysics scales with the length of the time step.

For example, a vertical velocity of 10 m/s for 10 s results in a 100‐m displacement of the air, resulting in a dry adiabatic cooling of around 1 K. Assuming that the air parcel starts at liquid saturation, this leads to SSLiq of 6.4% at +20 °C, increasing to 9.1% at −20 °C. SSIce also scales similarly with time step length. Because simulated convective updrafts can be as high as 50 m/s, this problem is not negligible even with rather short time steps.

These large SSLiq and SSIce values affect many processes within the microphysics, for example, the growth by vapor deposition of all ice hydrometeors scales linearly with SSIce. Additionally, the activation of CCN (using; Segal & Khain, ) occurs only where SSLiq ≥ 0, although the updraft velocity, rather than supersaturation, is used to determine which CCN are activated. The microphysical parameterization behaves very differently, and nonlinearly, depending on the input SSLiq and SSIce values. Repeatedly calculating microphysical process rates based on large supersaturations and then applying these rates over long time steps gives significantly different results compared to calculating the microphysical rates more frequently, and thus, with a smaller supersaturation.

By changing to parallel splitting, the results are much less sensitive to the model time step and the simulations are more consistent with what occurs in the atmosphere. With parallel splitting, all of the microphysical processes are calculated in an environment that is close to (i.e., never significantly exceeding) liquid supersaturation. In contrast, sequential‐update splitting always uses the SSLiq and SSIce values determined after the model dynamics. In convective updraft regions, and with long time steps, SSLiq could be several tens of percent. Consistently calculating microphysical process rates at such unrealistically high supersaturations is very inconsistent with the microphysical processes of the real atmosphere.

Sensitivity of Microphysical Process Rates to Model Time Step

The surface precipitation is substantially reduced in the default version of COSMO as the time step increases. However, the SUS+satad simulations and those with parallel splitting show much less time step dependence. To explain the differences, we performed a rain budget analysis. All processes that are a source or sink of rain water in the model were tracked and accumulated over the model simulation. This was done for each of the model simulations with different time steps, and for all three numerical splitting schemes in this paper.

The sources of rain mass are autoconversion of cloud water to rain, accretion of cloud water by rain, and melting of ice‐phase particles. The sinks of rain mass are freezing of rain to cloud ice, graupel, and hail; riming of rain on to any frozen particle; and evaporation of rain falling below cloud base in a subsaturated environment. Autoconversion and accretion are together labeled “warm rain” in Figure .

Rain water budget showing vertically integrated sources and sinks of rain water from the different physical processes represented in the model (a–c) and calculated contributions to surface precipitation based on these diagnostics (d–f). The eight columns for each process and calculation show the quantity for a simulations with a different time step, left‐to‐right: 1 to 15 s. The three rows show the process rates from simulations with sequential‐update splitting (top row), sequential‐update splitting with extra saturation adjustment within the model dynamics (middle row), and parallel splitting (bottom row).

Figures a–c show the vertically integrated sources and sinks for simulations with varying time steps and numerical splitting schemes. Within each category there are eight bars, representing different time step simulations (left‐to‐right: 1, 2, 3, 4, 5, 6, 10, 15 s respectively). Figures d–f show calculated quantities that reflect the relative importance of warm rain production, rain produced via the ice phase, the rain mass in the atmosphere at the end of the 2‐hr simulation, the diagnosed surface rain based on the rain water budget, and the actual precipitation output from the model. The precipitation total diagnosed from the rain water budget and actual precipitation agree to within 10%, and the trends with respect to changing model time step are very accurately reproduced. Detailed analysis of the geographical distribution of these processes (not shown) revealed that the rain water existing in the model at the end of the simulation has been recently created through autoconversion or accretion. Almost all of the evaporation occurs in regions directly underneath melting hydrometeors, and therefore, it is clear that almost all the surface precipitation occurs via ice‐phase processes. All snow and graupel melts before reaching the surface.

In additional to the vertically integrated process rates, the vertical profiles of the process rates are shown in Figure for simulations with 1, 5, and 15 s. These profiles are colored based on the process in Figure to which they contribute.

Vertical profiles of microphysical process rates contributing to the rain water budget in Figure . The lines are colored to coordinate with the process groups in Figure . The three rows show the process rates from simulations with sequential‐update splitting (a), sequential‐update splitting with extra saturation adjustment within the model dynamics (b), and parallel splitting (c). For each numerical setup, process rates for 1‐, 5‐, and 15‐s time step are shown. The right‐hand column shows the number concentration of ice particles produced through ice nucleation. For the 15‐s simulation in panel (a), the maximum value reaches 1.035 × 1011. The evaporation rates have been multiplied by 6 to show the differences more clearly.

When sequential‐update splitting is used (Figure , top row), there is a clear reduction of the rain produced through autoconversion and accretion (“warm rain”) as the time step is increased. This is the opposite of what would be expected from the simple model, where autoconversion rate is increased at longer time steps. The vertical profiles reveal that there is also a large change in the height at which the autoconversion and accretion occur (Figure a). Indeed, at lower heights, the autoconversion rate increases with increasing time step length as is expected from the simple model. However, there is a large reduction in both autoconversion and accretion higher in the cloud for longer time steps. The reason for this decrease is a very large increase in the freezing of cloud water (via ice nucleation; Figure a) at longer time steps, limiting the cloud water available for conversion to rain. As a result of the reduced source of rain water, the freezing of rain water and riming on to frozen particles is also reduced for the longer time steps (Figure a). The height at which riming of rain on to ice, graupel, and hail also changes substantially with time step (Figure a). There is no riming of rain on to snow because the two species never coexist in these simulations.

These changes can be traced back to a large sensitivity of the ice nucleation parameterization to the model time step. This impacts the other microphysical processes quite dramatically, partly because the ice‐phase processes (starting with ice nucleation) are calculated before the warm‐phase processes in the Seifert and Beheng () scheme. The ice nucleation parameterization used (Fletcher, ;  with modification from; Huffman & Vali, ; see  ) is explicitly dependent on the supersaturation with respect to ice and compared to the “maximum” supersaturation with respect to ice (that of water saturation, S₀). Therefore, with long time steps, ice nucleation becomes very active near the freezing level where $\frac{S_i}{S_0}>1$ can be achieved (or equivalently where SSLiq > 0). These numerous but small ice crystals quickly rime in the water‐rich environment producing a very large number of small graupel particles. These graupel particles are lofted high into the atmosphere by the strong updraft, but typically remain small. Therefore, these graupel particles do not contribute to precipitation, but instead, they remain within the atmosphere. As a result, the precipitation efficiency of the clouds is significantly reduced at longer time steps, and therefore, the surface precipitation is significantly reduced. At shorter time steps, the graupel particles remain near the updraft and collect supercooled water. Eventually, they grow large enough to be considered hail and continue to grow by collecting supercooled water. A large sensitivity to model time step of the riming of rain water on to hail (Figure a), due to reductions of both hail and rain water in the appropriate region, also contributes to the lack of precipitation reaching the surface.

In the lower two rows of Figure , the rain budget is shown for the SUS+satad simulations (Figure b) and those using parallel splitting (Figure c). The results for the rain budget are almost identical for the two setups, so they are described here together. In contrast to the default setup, there is an increase in the rain mass produced through warm rain processes for longer time steps. This increase in the warm rain process is almost exactly the same for the riming, and there is a slight decrease in the freezing of rain for the longest time steps. Despite the increase in riming, melting is not increased, which suggests that these extra rimed particles never fall to the melting level but instead stay in the atmosphere. Alternatively, the extra riming of rain water could be countered by a reduction of cloud riming, which was not tracked as part of the rain water budget. The evaporation of rain is approximately constant for all time steps. Therefore there is almost no sensitivity of the rain formed through the warm rain process alone (warm rain minus freezing and riming), and almost no change to rain from melting, as the model time step gets longer. Therefore, the total surface precipitation remains approximately constant across all model time steps. The only notable differences between the SUS+satad and parallel splitting process rates are the small increase in melting and the reduction of evaporation at long time steps when parallel splitting is used. These both lead to an increase in the melting‐minus‐evaporation term for the longest time steps when parallel splitting is used and a slight increasing trend of precipitation with increasing time step at the longest time steps.

Although the precipitation output from the SUS+satad and parallel split simulations are almost independent of time step, there remains a time step dependency in both the warm rain production and riming rates (Figures b and c, and b and c). In the SUS+satad simulations this is expected due to the time step dependence of the cloud water content produced by saturation adjustment, analogous to the simplified model in section . For the parallel split simulations, the time step dependence of the warm rain production is reduced by about half compared to the SUS+satad simulations (Figure c). This sensitivity remains as a result of the saturation adjustment that occurs after the microphysics during the previous time step. The time step dependence of the riming rates are explained by the changing availability of rain water. The reduction of the time step sensitivity of riming of rain on to hail under parallel splitting (Figure c), compared to sequential‐update splitting and SUS+satad is a key reason for the more consistent surface precipitation and hail totals for all model time steps.

The characteristics of the ice nucleation within the default COSMO setup helps explain why the sensitivities are much reduced in the SUS+satad or parallel splitting setups. Neither of these setups permit SSLiq to be diagnosed by the microphysics scheme. Therefore, the SSIce is limited, especially at higher temperatures and consequently the ice nucleation parameterization behaves consistently across time steps (i.e., there is no nucleation of ice at temperatures above −20 °C). This analysis of the rain budget does reveal that there are also other microphysical processes in the model that are also sensitive to the model time step. In simulations where the ice number concentration parameterization was changed to the Fletcher () parameterization (not shown), the time step sensitivity for surface rain was reduced from 40% to 14% for the sequential‐update simulations, the SUS+satad and parallel split simulations were unchanged.

Discussion

In this section, the broader importance of these findings is demonstrated. First, the impact of the choice of splitting method on the microphysical parameterization is discussed. Then, a short discussion of the impact of model time step when using the operational single‐moment microphysics scheme is presented. Finally, the impact of changing the operator splitting method on the robustness of scientific findings is explored.

Splitting Choices Within the Microphysics Parameterization

In addition to the splitting of the model dynamics and model physics, the individual processes within the microphysical parameterization are also implemented in succession using sequential‐update splitting. This approach ensures numerical stability; however, there are two drawbacks with this approach. First, like with the dynamics‐microphysics splitting, there is an inherent time step dependence by using the sequential‐update splitting approach. Second, the microphysical processes are calculated in a semiarbitrary order, meaning that any change to the order of the calculated processes will change the result of the microphysical parameterization calculation. Taking an extreme example, if all of the freezing processes occur at the end of the microphysical parameterization, then significant liquid water is likely to be present for riming processes; in contrast, with the freezing processes at the beginning (as in the Seifert & Beheng, , scheme), much of the liquid water could be frozen before the riming processes are calculated, which would severely reduce the amount of rimed ice produced in the model.

The dependence on both model time step and the order of the processes could be reduced by using parallel splitting for all process in the microphysics. However, to do so would make it possible that some water contents become negative because multiple processes could act to remove mass from a particular hydrometeor class. For example, if there are three separate processes, each which would remove 50% of the currently available liquid water, then applying parallel splitting would require 150% of the liquid water to be removed. Hence, the liquid water content would become negative and this numerical instability would result in the model crashing. Getting a numerically stable, positive‐definite, solution would require that each of the processes is clipped to maintain nonnegative mass contents for all the hydrometeor species. However, with so many separate processes involved this approach is complicated and also somewhat artificial.

Time Step Sensitivity With a Simpler, Single‐Moment Microphysics Scheme

To understand the importance of the time step sensitivity on simulations more typical of numerical weather prediction models, we replace the two‐moment Seifert‐Beheng microphysics scheme with the German Weather Service's operational single‐moment microphysics scheme (Baldauf et al., , section 2b). This microphysics scheme includes saturation adjustment before the microphysical processes are calculated. The simulations with this microphysics scheme for the same range of time steps is shown in Figure . Although the absolute precipitation total in the model domain is relatively constant for all the model time steps chosen (range of results within 15% of maximum value), the spatial distribution of precipitation is again strongly affected. For the shortest time steps, the heaviest precipitation (maximum: 20.5 mm) is in the center of the model domain and results from the first convective cell (initiated by the warm bubble). However, the longer time step simulations show much lower maximum precipitation from the initial convective cell (11.5 mm), and the largest precipitation totals are associated with the splitting cells later in the simulation and falls 25 km further east. Although the precipitation totals are quite similar, their magnitudes are higher and more concentrated on smaller regions when a shorter time step is used. These results are similar to that of the two‐moment scheme when saturation adjustment is applied during the model dynamics.

Similar to Figure , except for simulations using the single‐moment microphysics scheme from the operational COnsortium for Small‐scale MOdeling version. Only the precipitation is shown because hail does not exist in the one moment scheme. The y scale in the lower plot has increased.

Implication for Process Studies

The impact of model time step on the simulations is not only a matter of changing the precipitation amount or distribution—which could be corrected for by using postprocessing, for example. Additionally, the sensitivity of the model to changed forcings is also affected. As an example, it is clear from Figure , that both total precipitation amount and total hail amounts are sensitive to changes in the CCN concentration. However, it is also clear that the magnitude and sign of this sensitivity changes systematically as a function of model time step. The magnitude and sign of CCN effects on precipitation and hail also vary with model time step in the SUS+satad simulations (Figure ), although differently than in the sequential‐update splitting simulations. There is still a large uncertainty about the effect of aerosol concentrations on precipitating mixed‐phase cloud systems (Fan et al., ). A dependence of the model results to the model time step could be contributing to this uncertainty.

Additionally, we have tested the sensitivity of total surface hail to changing ice nucleation rate. In these simulations the default ice nucleation rate was scaled up or down by up to a factor of 100, and simulations were run in three different configurations: (a) default COSMO setup with sequential‐update splitting and a 2‐s time step, (b) with SUS+satad with a 6‐s time step, and (c) with parallel splitting and a 6‐s time step.

As shown in Figure , the absolute value and sensitivity of the hail fall to the ice nucleation rate is markedly different for each of the three model setups. In the setup with parallel splitting, the sensitivity to the ice nucleation scaling is the largest, while the sensitivity is smaller and nonmonotonic for the other two model setups. Thus, depending on the setup of the model that is being used, vastly different scientific conclusions could be drawn about the importance of the ice nuclei concentration for hail precipitation.

Total hail fall reaching the surface during the first 6 hr from simulations where the ice nucleation rate is scaled from 0.01 to 100. Each line is for a different setup of the model numerics. The simulations marked by the orange line included saturation adjustment within the model dynamics. All model setups included saturation adjustment after the microphysics.

Conclusion

In this paper, we have shown that a weather prediction model coupled to a two‐moment microphysics scheme produces significantly different results when the model time step is varied between 1 and 15 s. Simulations with the shortest time steps have significantly more precipitation, a different location for the heaviest precipitation and a different sign for the affect of changing CCN on precipitation than the longest time step simulations. Such a time step sensitivity is especially important for numerical weather prediction simulations, which use long time steps to decrease computational cost.

The sensitivity arises because of the use of sequential‐update splitting between the model physics and dynamics. This is particularly important in convective updraft areas, where large SSLiq is created and subsequently used to drive the microphysical processes. However, when a hybrid parallel splitting is used or an additional saturation adjustment is added inside the model dynamics, the time step sensitivity is reduced, because the large SSLiq values are no longer passed directly to the microphysics. In the sequential‐update splitting setup with the longest time step, this supersaturation results in a large increase of ice nucleation rate at relatively high temperatures. Due to an increase in the number of ice particles and a reduction in the available supercooled water, the convective clouds in these simulations are not as efficient at converting condensed water to precipitation, for which riming is an important process. Consequently, the precipitation and hail is decreased as the time step is increased.

The reductions of precipitation and hail totals with increasing time step are less severe in both modified setups than the default setup. However, neither of the additional setups provide fully time step independent simulations either. The intensity and location of the largest accumulated precipitation also reduce with increasing time step in the SUS+satad setup, but they are maintained for all time steps when parallel splitting is used. The changing sign of CCN sensitivity also exists for the SUS+satad simulations, but when using parallel splitting the impact of increasing CCN is to systematically reduce the precipitation and hail totals for all time steps used. The process rate analysis reveals that the autoconversion and accretion rates still change with model time step under SUS+satad and parallel splitting, as do the riming rates. The causes for these differences needs to be further understood in future work.

Our findings that parallel splitting is less sensitive to the model time step than sequential‐update splitting contradicts several published papers which show advantages of sequential‐update splitting (e.g., Beljaars et al., ; Dubal et al., ; Rood, ; Williamson, ). The reason for this is probably that there is significant vertical velocity in our simulations, which produces significant SSLiq and significantly affects the microphysics. In the other studies, there were no sink terms in the model that explicitly depend on the input model state and hence the time step length. Therefore, the microphysics of convection‐permitting simulations needs to be treated differently.

There are numerous negative consequences of this time step dependence. First, and most obvious from this study, is the potential for contradictory findings from process studies when choosing a different model time step. Taking the example of the effect of CCN on precipitation and hail, all effects from a large increase to a large decrease are possible depending on the chosen time step and process splitting methodology.

This uncertainty regarding the impact of CCN on convection is in agreement with numerous published, contradictory results about the impact of aerosol on clouds. For example, two studies for the same day found opposite conclusions about the effect of CCN (Khain et al., ; Noppel et al., ). While these results could be due to the different models used in these two studies, the effect of model sensitivity to time step could be affecting one or both of the studies.

The authors do not have access to the numerical implementation of all models, but many models are typically built in a similar way and are therefore likely to have similar problems. We have run simulations using the WRF model including the two‐moment NSSL microphysics scheme and found that there are also significant changes of the simulated storms when the model time step changes. However, the WRF simulations changed less systematically with the model time step than those of the COSMO model. This supports the discussion in Lebo et al. () that “simulated results can differ both quantitatively and qualitatively when switching from a fixed timestep to an adaptive timestep method” in the WRF model. Additionally, different time integration used in ECHAM‐HAM climate model used to solve the sulfuric acid gas evolution equation led to substantially different concentrations and aerosol nucleation rates (Wan et al., ). Recent improvements in the ECMWF model came from changing the implementation of microphysical processes to make the model more robust to time step (Forbes, ). The extent to which other models and other physical processes are affected by similar issues is unknown. Nevertheless, the existence of a strong effect of the model time step in these state‐of‐the‐art models suggests that it is likely an issue that affects other models and deserves further investigation.

Acknowledgments

The research leading to these results has been done within the subproject “Microphysical Uncertainties in Deep Convective Clouds and their Implications for Data Assimilation” of the Transregional Collaborative Research Center SFB/TRR 165 “Waves to Weather” (www.wavestoweather.de) funded by the German Research Foundation (“Deutsche Forschungsgemeinschaft”, DFG). This work was performed on the computational resource ForHLR II funded by the Ministry of Science, Research and the Arts Baden‐Württemberg and DFG. Data from the model simulations are available through KIT's open‐data portal (doi:10.5445/IR/1000084513). We thank the reviewers for their contributions to significantly improving this manuscript.