Geoscientific Model Development Analyzing numerics of bulk microphysics schemes in community models : warm rain processes

Implementation of bulk cloud microphysics (BLK) parameterizations in atmospheric models of different scales has gained momentum in the last two decades. Utilization of these parameterizations in cloud-resolving models when timesteps used for the host model integration are a few seconds or less is justified from the point of view of cloud physics. However, mechanistic extrapolation of the applicability of BLK schemes to the regional or global scales and the utilization of timesteps of hundreds up to thousands of seconds affect both physics and numerics. We focus on the mathematical aspects of BLK schemes, such as stability and positive-definiteness. We provide a strict mathematical definition for the problem of warm rain formation. We also derive a general analytical condition (SMcriterion) that remains valid regardless of parameterizations for warm rain processes in an explicit Eulerian time integration framework used to advanced finite-difference equations, which govern warm rain formation processes in microphysics packages in the Community Atmosphere Model and the Weather Research and Forecasting model. The SMcriterion allows for the existence of a unique positive-definite stable mass-conserving numerical solution, imposes an additional constraint on the timestep permitted due to the microphysics (like the Courant-Friedrichs-Lewy condition for the advection equation), and prohibits use of any additional assumptions not included in the strict mathematical definition of the problem under consideration. By analyzing the numerics of warm rain processes in source codes of BLK schemes implemented in community models we provide general guidelines regarding the appropriate choice of time steps in these models.


Introduction
In the last decade there has been only one study (Morrison and Gettelman, 2008) (MG08) that discusses explicit time integration scheme used to advance governing microphysical equations in a bulk microphysics (BLK) scheme.MG08 implemented a two-moment BLK scheme with a diagnostic treatment of precipitating hydrometeors in the Community Atmosphere Model (CAM) version 3 (Collins et al., 2006).To advance governing microphysical prognostic equations in MG08, time splitting is applied to separate sedimentation from the rest of the microphysical processes, and rain and snow sedimentation are treated diagnostically similarly to Ghan and Easter (1992).For cloud water and cloud ice sedimentation, an upstream advection scheme and an explicit time integration scheme are used.The Courant-Friedrichs-Lewy condition is always satisfied because the microphysical time step is sub-divided into smaller equal time steps.This procedure assures stability and positive definiteness.To advance other microphysical equations, an explicit Eulerian time integration framework is used.Because stability and positive-definiteness criteria are not defined, sometimes hydrometeors' mixing ratio might be negative.To prevent negative mixing ratios, a so-called "mass conservation" technique is used to calculate artificially modified growth rates for some microphysical processes but the time step remains unchanged.We speculate that the "mass conservation" approach can influence such cloud characteristics as effective radius, radar reflectivity, precipitation fluxes, and cloud radiative properties as well as the amount of accumulated precipitation and its spatial and temporal distributions.We also speculate that when applied in global scale earth system models with time steps of 10-20 min, a non-stable Published by Copernicus Publications on behalf of the European Geosciences Union.
and non-positive-definite time integration scheme could influence the global water cycle causing artificial precipitation patterns that are then used as surface boundary conditions for ocean, land, lake, and sea ice models.
In the Weather Research and Forecasting (WRF) model version 3 (Skamarock et al., 2008), both the single moment BLK schemes and the Morrison-Curry-Khvorostyanov double-moment BLK scheme (Morrison et al., 2005) (MCK05) might share similar deficiencies of non-positive and unstable solutions for warm rain processes if the microphysical time step used is greater than a few tens of seconds.This feature of BLK schemes implemented in community models (CAM and WRF) could lead to possible erroneous conclusions regarding the role of cloud microphysics and their influence on radiation or dynamics (amongst others) when relatively long time steps are used for integration.
To avoid the uncertain performance of BLK microphysics schemes, if relatively long time steps are used for model integration and to improve the creditability of precipitation amount calculations, we derive the necessary condition (referred to as SM-criterion) to keep an explicit Eulerian time integration scheme stable and positive-definite regardless of the parameterizations used for warm rain processes in BLK schemes.The SM-criterion constitutes the existence of a unique positive-definite stable numerical solution and imposes constraints on the time step permitted (like the Courant-Friedrichs-Lewy condition for the advection equation).We highlight that in addition to the limitations on the time step imposed by the dynamics, there also exists a limitation due to the microphysics.We also define well-behaved, conditionally well-behaved, and non-well-behaved Explicit Eulerian Bulk Microphysics Code (EEBMPC) classes and show that source codes of BLK schemes, which were originally developed for use in cloud-resolving models and implemented in community models, belong to conditionally wellbehaved EEBMPC class.We also provide recommendations regarding integration time steps for prospective simulations with WRF and CAM.
The paper is organized as follows.The general considerations are given in Sect. 2. The growth rate calculation due to warm rain processes in BLK schemes are discussed in Sect.3. The comprehensive multi-step analysis of numerics for the system of differential equations that governs processes of warm rain formation in BLK scheme is presented in Sect. 4. The first step in our analysis is the derivation of a mass-conserving positive-definite analytical solution for the linearized differential-difference equations, presented in Sect.4.1.The second step in our analysis is the derivation of a numerically explicit Eulerian solution for the finite-difference equations, outlined in Sect.4.2.The third step in the comprehensive analysis is the stability analysis for the numerical explicit Eulerian solution, given in Sect.4.3.In Sect.4.4 we show how the utilization of the "mass conservation" technique might cause the violation of the stability and positive-definiteness conditions when analytical representation of autoconversion and accretion growth rates are known, using one of the BLK schemes as example.Discussion and recommendations are provided in Sect. 5.

General consideration
We consider the following system of equations for bounded positive X(t) and Y (t) with initial conditions where F (X) and G(X, Y ) are both positive and bounded, and positive τ is given.
We are looking for a numerical solution for X(n + 1) and Y (n + 1) with initial conditions X(n) = X 0 and Y (n) = Y 0 such that at each time step "n" satisfies the conservation equation as well as positiveness and boundedness conditions where X(n) and Y (n) are X and Y at the beginning of time step "n".Explicit finite-difference scheme with time step τ is written as It is clearly seen that this solution conserves the sum of X and Y .By adding expressions ( 6)-( 7) we get finite-difference analog for the conservation equation given by Eq. (3): The values of X(n + 1) and Y (n + 1) at the next time step are always bounded and positive only if Expression (8) determines the time step permitted to keep the solution ( 6)-( 7) bounded and positive.To solve the system we have to specify X, Y , F (X), G(X, Y ), and τ as well as determine a stability condition.

Warm rain processes in community models
Applying these general considerations to cloud physics and defining we obtain the following system of equations that governs the process of "warm" rain formation in prognostic BLK schemes: where Q c is cloud water mixing ratio, Q r is rain water mixing ratio, and PAUTO and PACCR are their changes due to auto-conversion and accretion, respectively.By adding Eqs. ( 9) and ( 10), we get Equation ( 11) has the simple physical meaning.For the "warm" rain formation process total water mixing ratio Different single-moment and double-moment BLK schemes used in CAM and WRF and referenced thereafter as RaschCAM (Rasch and Kristjansson, 1998) and KESSLER (Kessler, 1969), LIN (Lin et al., 1983), ETA (Ferrier, 1994), TAO (Tao et al., 2003), THOMPSON (Thompson et al., 2004), MORRISON (Morrison et al., 2005), and WSM6 (Hong and Lim, 2006) schemes, respectively, calculate auto-conversion and accretion growth rates using a different non-linear analytical representation for functions PAUTO and PACCR.The system of non-linear differential Eqs. ( 9)-( 10) is linearized and solved using an explicit Eulerian (EE) time integration scheme.Our analysis of source codes of MG08 scheme in CAM and TAO, THOMPSON, MORRISON, and WSM6 schemes in WRF, respectively, reveals that positiveness criterion (SM-criterion) similar to the inequality (8) and given by is never checked.
In our analysis of source codes for the BLK schemes referred to above, we use the following definitions.Based on the validation of the SM-criterion we introduce definitions for three different classes that are a well-behaved EEBMPC class, a conditionally well-behaved EEBMPC class, and poorly-behaved EEBMPC class.
If the SM-criterion has always been validated and respected in EEBMPC, we define such a code as belonging to well-behaved EEBMPC class.The remarkable feature of a well-behaved EEBMPC is that it assures a unique stable positive-definite mass-conserving numerical solution (referred to as correct solution thereafter) for the governing differential equations.To the best of our knowledge, there is no well-behaved EEBMPC implemented in Community models.
If the SM-criterion has never been checked in EEBMPC, we define such a code as belonging to a conditionally well-behaved EEBMPC class or a poorly-behaved EEBMPC class, whose common feature is that both rely on a so-called "mass conservation" technique in an attempt to avoid negativeness of hydrometeors' mixing ratios and make an explicit Eulerian time integration scheme positive-definite.As it is shown in the next section, the quintessence of the "mass conservation" technique is an incorrect assumption: it assumes the existence of a positive numerical solution for the governing differential equations in a time interval where this solution is not defined.
The distinguishable feature of a conditionally wellbehaved EEBMPC is that the SM-criterion is satisfied even if it has never been checked.As used in cloud resolving models or large eddy simulation models with temporal resolution about a few seconds, conditionally well-behaved EEBMPC provides a correct solution for governing differential equations because the limitation on the time step imposed by dynamics is more restrictive than that imposed by microphysics (SM-criterion), and, in fact, the "mass conservation" technique is never applied.The transition between conditionally well-behaved EEBMPC and poorly-behaved EEBMPC is determined by the SM-criterion.If the time step used in a host model (CAM or WRF) is as those typically used for regional scale and especially large scale simulations and the SM-criterion is occasionally violated, a conditionally well-behaved EEBMPC would become a poorlybehaved EEBMPC.The eventual feature of a poorly-behaved EEBMPC is that it does not provide a correct solution for governing differential equations.
To avoid the use of an artificial "mass conservation" technique, mimic the utilization of well-behaved EEBMPC, and provide guidelines regarding time steps permitted for model integration for prospective CAM and WRF users, we present maximal time steps for autoconversion and accretion processes as well as effective maximal time steps for both processes that are necessary to keep stability and positivedefiniteness of an explicit Eulerian time integration scheme.These maximal time steps, which satisfy SM-criterion, are calculated for specified values of cloud water and rain water mixing ratios and two values of cloud droplet concentrations (N c ) for different BLK schemes and are shown in Table 1.The rightmost column in Table 1 shows the maximal time step permitted, if a particular BLK scheme is chosen, for a simulation.The values in paretheness correspond to different Table 1.Maximal time steps permitted to keep explicit Eulerian time integration scheme stable and positive-definite for autoconversion, accretion, and due to both processes in CAM and WRF bulk microphysics schemes.For comparison, the same values for Beheng (1994), Seifert and Beheng (2001), and Seifert and Beheng (2006) parameterizations are presented in gray.Maximal time steps shown for Q cloud droplet concentrations (N c ) and indicate dependence of autoconversion or accretion on N c = 10(100) cm −3 .For example, the KESSLER scheme representation of both autoconversion and accretion do not depend on N c , and for the KESSLER row, values of time step are equal in each column.
Non-dependence on N c , in the autoconversion column, manifests non-applicability of KESSLER, LIN, and TAO schemes for cloud-aerosol interaction simulations, whereas, only the THOMPSON scheme accounts for dependence of accretion on N c .Except for this scheme, all other WRF BLK schemes under consideration can be used for regional scale simulations if the time step in the host model does not exceed two to three hundred seconds.However, it should be noted that this conclusion is valid only for specific Q c = 1.0 g kg −1 and Q r = 0.5 g kg −1 used to calculate maximal time steps presented in Table 1.
To demonstrate the instantaneous dependence of the effective maximal time step on typical cloud water mixing ratio and rain water mixing ratio for different cloud types we define a SM-number (N sm ) as It should be noted that there is an obvious relationship between the SM-criterion and the SM-number.The SM-criterion is valid (an explicit Eulerian finite-difference scheme is stable and positive-definite) if Thus, the maximal time step permitted to keep an explicit Eulerian time integration scheme stable and positive-definite corresponds to Maximal time steps calculated according to expression (15) for TAO, THOMPSON, MORRISON, and WSM6 WRF BLK schemes as functions of Q c and Q r for two different droplet concentration N c = 10 cm −3 and N c = 100 cm −3 , which are used as a proxy for clean "maritime" and "continental" clouds, are shown on Figs. 1 and 2 and Figs. 3 and 4, respectively.For "clean maritime" clouds, Fig. 1 shows instantaneous dependence of maximal time step on Q c for Q r = 0.1 g kg −1 (top left), 0.5 g kg −1 (top right), 1.0 g kg −1 (bottom left), and 3.0 g kg −1 (bottom right), respectively, whereas Fig. 2 shows instantaneous dependence of maximal time step on Q r for Q c = 0.1 g kg −1 (top left), 0.5 g kg −1 (top right), 1.0 g kg −1 (bottom left), and 3.0 g kg −1 (bottom right).Figures 3 and 4 replicate Figs. 1  and 2, respectively, for "clean continental" clouds.
The set of these four figures represents a simple yet powerful tool to analyze the behavior of a BLK microphysics scheme.Since a conventional way to validate new microphysics parameterizations is through the use of a single column model (SCM), observations (vertical profiles of cloud water mixing ratio, rain mixing ratio, and cloud droplet concentration) and data from Figs. 1-4 can be used to analyze theoretical vertical profiles of SM-criterion.These vertical profiles provide a useful way to make an appropriate choice for the time step used for a SCM integration instead of arbitrary values as is conventionally done.Moreover, in the case of 2-D or 3-D simulations, the vertical profiles of the SMcriterion show additional limitations imposed on the time step used in a multidimensional host model.For example, limitation on the time step provided in the WRF User Guide is given as where τ wrf is the time step in seconds and X wrf is the spatial resolution in kilometers.However, this constraint does not account for additional restrictions due to the microphysics imposed by the SM-profiles that ensure an explicit Eulerian scheme stability and positive definiteness and reliability of model output.In fact, for regional or large scale WRF simulations with a time step chosen according to inequality ( 16) violation of the SM-criterion at different times, altitudes, and spatial locations leads to a non-positive-definite numerical solution for the governing warm rain differential equations.
To show that the SM-criterion is a necessary condition for an explicit Eulerian time integration scheme stability and positive-definiteness and to demonstrate the consequences of using the "mass conservation" technique, the next section describes the multi-step analysis used.

Analytical and numerical solutions and stability analysis
Different single-moment and double-moment BLK schemes used in community models formulate auto-conversion PAUTO and accretion PACCR growth rates in a variety of ways providing different non-linear functional dependences on Q c , Q r , and N c .The system of nonlinear differential Eqs. ( 9)-( 10) that governs the processes of warm rain formation can be solved only numerically using iterative methods.However, if some linearization is assumed, they could also be solved analytically.
The analytical solution ( 23)-( 24) is bounded and positivedefinite if and only if These inequalities determine the maximal time step permitted to keep a bounded and positive solution.Both are satisfied if the SM-criterion is valid: Thus, the SM-criterion determines the sufficient and necessary positiveness condition for the analytical solution to the system of differential Eqs. ( 21)-( 22) regardless of the specific formulations for autoconversion PAUTO and accretion PACCR growth rates.Condition ( 27) also determines the applicability of the linearization given by expressions ( 17)- (20).At this point it should be clear that any assumption regarding the existence of the analytical solution for a time greater than that given by condition ( 27) is not sensible mathematically.The analytical solution for the linearized differential-difference Eqs. ( 21)-( 22) permanently exists for any "t" on time interval 0 ≤ t ≤ τ max only.

Warm rain processes: explicit Eulerian time integration scheme
The finite-difference analog for the system of nonlinear differential Eqs. ( 9)-( 10) that govern processes of warm rain formation can be given as where q n c and q n+1 c , q n r and q n+1 r are initial and new values of cloud water mixing ratio and rain mixing ratio, respectively.Time representations for auto-conversion and accretion growth rates are still not specified.In general, Eqs. ( 28)-( 29) can be solved only using iterative numerical methods that need significant computational time.However, if some linearization is assumed, non-iterative computationally efficient numerical methods can be used.For example, in an explicit Eulerian scheme a linearized explicit form is used for both auto-conversion and accretion growth rates: www.geosci-model-dev.net/5/975/2012/Geosci.Model Dev., 5, 975-987, 2012 where C n u and C n a are given by Explicit representation means that both auto-conversion and accretion growth rates can be calculated at the beginning of the microphysical time step because q n c and q n r are known.With expressions (30)-(33), Eqs. ( 28)-( 29) are as follows: Solving for q c and q r we get a numerical solution for linearized differential Eqs. ( 9)-( 10): It is clearly seen that this solution conserves mass.By adding expressions ( 36)-( 37) we get the finite-difference analog for the mass conservation equation given by Eq. ( 11): Despite the fact that the solution ( 36)-( 37) conserves mass, it is not positive-definite.Whereas q n+1 r is always positive, q n+1 c sometimes might be negative.The numerical solution (36)-( 37) is bounded and positive-definite if and only if These inequalities determine maximal time step permitted to keep a bounded and positive numerical solution.Both inequalities are satisfied if the SM-criterion is valid: Thus, the SM-criterion provides the necessary condition for the explicit Eulerian finite-difference scheme ( 34)-( 35) to be positive-definite regardless of the parameterization formulae used for autoconversion and accretion growth rates ( 30)-( 33).
An observation that the solution ( 23)-( 24) for differentialdifference equations and the solution (36)-(37) for finitedifference equations coincide is important, and its mathematical meaning is that the finite-difference scheme is stable for fixed timesteps that do not exceed the maximal timestep given by the SM-criterion for the finite-difference equations (40).It should be clear that any attempt to solve the finite-difference Eqs. ( 34)-( 35) using a timestep that is greater than that given by ( 40) is neither mathematically nor physically sensible because this situation is not governed by these equations in an explicit Eulerian time integration framework.For a different time integration framework, the positive-definiteness condition and stability condition might differ.Stability is a very important issue that makes the finitedifference equations different from the differential-difference equations for which the stability problem is not relevant.The analysis of stability is of crucial importance for any finitedifference scheme and should be done before its implementation in a numerical model.

Stability analysis: explicit Eulerian scheme
The numerical solution ( 36)-( 37) can be written as the matrix equation: The matrix characteristic equation for system (41) has the following form: For the finite-difference scheme to be stable it is necessary that all roots λ 1,2 of its characteristic equation satisfy However, in the case the numerical solution (36)-(37) might oscillate.This fact contradicts the conditions of positiveness ( 38)-(39).Thus, instead of inequality (43) λ 1,2 must satisfy To find stability conditions for the scheme given by expressions (32)-( 35) all the roots of the matrix characteristic equation (42) have to be found.We then get the algebraic characteristic equation The solutions are as follows The first root λ 1 is given by and stability inequality (45) for λ 1 is always satisfied.The second root λ 2 is given by According to (45) it must satisfy the following inequality The right inequality is held unconditionally, but for the left inequality to be valid it is necessary that The condition given by expression ( 51) is necessary for the computational stability of the finite-difference scheme given by expressions (32)-( 35).Observation that conditions (40) and ( 51) coincide permits us to conclude that the SMcriterion provides the necessary condition for the explicit Eulerian finite-difference scheme given by Eqs. ( 34)-( 35) to be stable and positive-definite regardless of the parameterization used for autoconversion and accretion growth rates.
Because the stability of an explicit Eulerian time integration scheme for microphysical governing equations used in BLK schemes has never been discussed, the effect of the violation of the SM-criterion on its stability and positive definiteness is hidden.Since the validation of the SM-criterion is not reproduced in EEBMPCs used in community models, these codes belong to the conditionally well-behaved EEBMPC class.Thus, we conclude that if relatively long time steps are used for WRF integration and the SM-criterion is occasionally violated, the source code implementations of the TAO, THOMPSON, MORRISON, and WSM6 schemes in the official WRF distribution would belong to poorlybehaved EEBMPC that do not provide a correct numerical solution for governing differential equations.
Although the source codes for these four schemes share the same deficiencies, in the following section we provide a detailed analysis of numerics of warm rain processes in conditionally well-behaved EEBMPC for the MORRISON scheme in WRF (Morrison et al., 2005) as well as Morrison-Gettelman scheme in CAM (Morrison and Gettelman, 2008), because the numerical treatment of cloud water mixing ratio in both schemes is identical.

Warm rain processes in WRF: MCK05 numerical solution
The finite-difference analog for the system of differential equations for the double-moment BLK scheme is not presented and discussed in Morrison et al. (2005).
In this scheme the warm rain formation processes are governed by the system of differential equations (Khairoutdinov and Kogan, 2000) (KK2000): Using "reverse engineering" (translating source code documented in WRF to scientific notation used in the theory of finite-difference schemes), the finite-difference analog for the system ( 52)-( 53) can be written as where explicit representations for both PAUTO and PACCR are used: Explicit representation means that both the auto-conversion PAUTO and accretion PACCR growth rates, which are calculated at the beginning of the microphysical time step using known q n c , q n r , and N n c , are constants and can not be "adjusted".
Solving for q n+1 c and q n+1 r we get a numerical solution for the differential Eqs. ( 52)-( 53): It is clearly seen that this solution conserves mass.By adding expressions ( 58)-( 59) we get a finite-difference analog for the mass conservation equation given by Eq. ( 11): Although the solution ( 58)-( 59) conserves mass, it is not positive-definite, whereas q n+1 r is always positive, q n+1 c sometimes might be negative because the positiveness condition given by the SM-criterion for the MORRISON scheme where τ max is the time step permitted to keep positive solution, is not satisfied.
To avoid negative q n+1 c , similar to the approach usually employed in other BLK schemes (e.g., Reisner et al., 1998) reduced artificial auto-conversion AAUTO and accretion AACCR rates are used (through the "mass conservation" technique): The expressions for the "adjusted" growth rates ( 61) and ( 62) deserve some comments.First, their analytical representation is quite different from those provided by KK2000.Second, their explicit dependence on timestep τ not found in the www.geosci-model-dev.net/5/975/2012/Geosci.Model Dev., 5, 975-987, 2012 original formulae ( 56)-( 57) is clearly seen.Third, it is assumed that these expressions remain valid during timestep τ provided by a host model.However, as explained above, any attempt to use a timestep that is greater than that given by the SM-criterion ( 60) is not relevant in an explicit Eulerian time integration framework.
Because the SM-criterion is never checked in the EEBMPC for the MORRISON scheme, we classify this code as belonging to the conditionally well-behaved EEBMPC class.However, if the SM-criterion is violated for a particular "set" of {q c ,q r ,N c , and τ } passed by a host model, the source code for the MORRISON scheme would become poorlybehaved EEBMPC.An attempt to avoid negativeness of q n+1 c calculated according to the finite-difference Eq. ( 58) by applying a "reduced" autoconversion AAUTO and accretion AACRR given by ( 61) and ( 62), respectively, that act during timestep τ > τ max (the "mass conservation" technique) is artificial and has nothing in common with the numerical solution for the differential Eqs. ( 54)-( 55) using the explicit Eulerian finite-difference Eqs. ( 58)-( 59) that have no positivedefinite and stable solution for τ > τ max .
If the SM-criterion is not respected, a poorly-behaved EEBMPC in the MORRISON scheme creates a virtual microphysics reality characterized by a "virtual" cloud water mixing ratio (q n cv ) and rain water mixing ratio (q n rv ) that are used instead of a "real" cloud water mixing ratio (q n c ) and rain water mixing ratio (q n r ) supplied by the host model.These virtual numbers can easily be calculated using the following procedure.If for input {q n c ,q n r ,N n c , and τ } supplied by a host model, N sm > 1, artificially "adjusted" AAUTO and AACRR rates are calculated using formulae (61) and ( 62).Then a system of two equations for "virtual" q n cv and q n rv derived by a) substitution of q n cv and q n rv instead of "real" q n c and q n r , respectively, in Eqs. ( 56)-( 57) and b) replacement of PAUTO and PACRR with AAUTO and AACRR in (56) and ( 57), respectively, has to be solved: The remarkable feature of these "virtual" solutions q n cv and q n rv is that a "virtual" SM-number for the MORRISON scheme (N mv ) defined as is always equal to one.
The mathematical meaning of the equality is that there is more than one available way to "adjust" original growth rates.The only condition is that a sum of artificially "adjusted" positive virtual autoconversion VAUTO and accretion VACCR growth rates has to be equal to q n c /τ .Thus, instead of the equality (66) a general definition for the "virtual" SM-number (N smv ) that remains valid regardless of the parameterization for autoconversion and accretion processes can be written as However, even if a) virtual "adjusted" VAUTO and VACCR rates can be calculated by different methods (for example, for the MORRISON scheme these rates are given by AAUTO and AACCR, respectively) and b) their sum can exactly be equal to q n c /τ , there is no additional equation similar to (58) that can be used for the calculation of q n+1 c .As our analysis in Sects.4.1-4.3shows, the finite-difference Eq. ( 58) is not valid for an arbitrary chosen timestep τ > τ max or, equivalently, N sm > 1.
If N sm > 1, the expression (67) would constitute an additional constraint that is not included in the strict mathematical definition of the problem introduced in Sect. 2. This expression is the general "mathematical foundation" of the "massConservation" technique, whose quintessence is an incorrect assumption: it assumes the existence of a positivedefinite explicit Eulerian numerical solution in a time interval that is greater than that given by the general condition N sm ≤ 1, whereas our analysis shows that such a solution does not exist in an interval where τ > τ max .
The physical meaning of the general virtual SM-number (N smv ) given by ( 67) is that "real" cloud water is completely depleted (q n+1 c = 0) by "adjusted" rates acting during timestep τ (according to the finite-difference equation ( 58)), and it is thought, that the problem of negative cloud water mixing ratio on the next timestep is eliminated.The artificial growth rates, that use "virtual" q n cv and q n rv (for example, for the MORRISON scheme these virtual numbers are given by the solution for the system (63)-( 64)), are then passed to a host model for post-processing analysis.
Although our analysis is focused on warm rain processes, we highlight that the inclusion of an ice phase makes the SM-criterion more restrictive because additional solid hydrometeors compete for the available cloud water.For example, the source code implemented into CAM (Morrison and Gettelman, 2008) and GFDL AM3 GCM (Salzmann et al., 2010) uses a diagnostic treatment of precipitating hydrometeors.However, the numerical treatment of cloud water generalized to include ice-phase remains similar to the one discussed above.To avoid numerical instability two equal time substeps are used, and to avoid non-positivedefiniteness "mass conservation" is applied to modify original growth rates due to numerous warm, ice, and mixed microphysical processes.Our comprehensive analysis relevant to warm processes shows that these growth rates are affected if a more restrictive criterion than the SM-criterion is not satisfied for timesteps typically used in GCM simulations (10-20 min).Then, the transition from conditionally wellbehaved EEBMPC to poorly-behaved EEBMPC occurs.In this case, the output of a poorly-behaved EEBMPC contains artificially modified growth rates due to different microphysical processes.

Discussion
To date, there are no studies that analyzed the numerics of bulk microphysics schemes with prognostic treatment of precipitating hydrometeors implemented in WRF (TAO, THOMPSON, MORRISON, and WSM6).Moreover, a finite-difference analog for each of these schemes has never been provided.Our analysis of source codes for these schemes in WRF reveals that a non-positive-definite explicit Eulerian time integration scheme is used to advance finitedifference microphysical equations.
Focusing on the mathematical aspects of BLK schemes, such as stability and positive-definiteness, we provide a strict mathematical definition for the problem of warm rain formation.We derive a general analytical condition (the SMcriterion) that remains valid regardless of parameterizations for autoconversion and accretion processes in an explicit Eulerian time integration framework used to advanced finitedifference equations that govern warm rain formation processes in BLK microphysics schemes.We also prove that the SM-criterion is a necessary condition of positive definiteness for the analytical solution of the linearized equations as well as a necessary condition of stability and positive-definiteness for an explicit Eulerian time integration scheme.The SMcriterion constitutes the existence of a unique positivedefinite stable mass-conserving numerical solution, imposes an additional constraint on the timestep permitted due to the microphysics (like the Courant-Friedrichs-Lewy condition for the advection equation), and prohibits the use of any additional assumptions not included in the strict mathematical definition of the problem under consideration.In general, preciseness of numerical scheme and time truncation errors would also be numerical issues to consider if there is a proof that a numerical scheme is stable and positive-definite.
Our analysis in Sects.4.1-4.3shows the non-existence of a unique positive-definite stable numerical solution in an explicit Eulerian time integration framework for differential equations that govern warm rain microphysical processes for microphysical environmental conditions and an arbitrary chosen timestep for which N sm > 1.This conclusion implies that any additional assumptions not included in the strict mathematical definition of the problem are not valid.One of these additional assumptions is the extrapolation of the existence of a positive-definite explicit Eulerian numerical solution to a time interval that is greater than that given by the general SM-criterion.The latter assumption is the quintessence of the so-called "massConservation" technique that assumes that "adjusted" growth rates are applicable in a time interval where a positive-definite numerical solution does not exist.An understanding of the fact that the numerical solution does not exist on an arbitrary chosen time interval is sufficient to reject the utilization of the "mass conservation" technique, and any additional proof for its rejection is not needed.
We highlight that the utilization of "mass conservation" technique applied to warm rain processes in an explicit Eulerian time integration framework is an incorrect attempt to avoid negativeness of cloud water mixing ratio.The "mass conservation" approach is conceptually incorrect because it relies on an assumption that "reduced" with respect to "original" autoconversion and accretion growth rates act during a given time step.This assumption contradicts a general rule used for the derivation of an explicit Eulerian finite-different representation for governing differential equations.In an explicit Eulerian framework, the "original" growth rates are known constants calculated at the beginning of each microphysical time step and can not be changed.
It is conventionally thought that the WRF model can be applied to a broad range of spatial scales from large eddy up to global scale simulations.However, for a prospective WRF simulation at a regional or global scale, the time step chosen according to recommendation provided in the user guide can cause occasional violations of the SM-criterion at different times, altitudes, and spatial locations.An inappropriate choice of the time step leads to non-positive-definite numerical solution for the microphysical governing differential equations and degradation of the ability of WRF to calculate precipitation amount and its spatial and temporal distribution.By introducing the concept of the SM-criterion vertical profile we provide a simple yet powerful tool that permits a rough estimation of an additional limitation imposed on the time step by microphysics and an appropriate choice for a time step for WRF simulations.
Depending on the validation of the SM-criterion in an EEBMPC, we introduce a definition for wellbehaved EEBMPC, conditionally well-behaved EEBMPC, and poorly-behaved EEBMPC.In a well-behaved EEBMPC, the SM-criterion is always validated and satisfied, and a remarkable feature of well-behaved EEBMPC is an assurance of correctness of a numerical solution for the governing differential equations.If the SM-criterion is never validated, EEBMPC is assigned to a conditionally well-behaved EEBMPC class or poorly-behaved EEBMPC class, whose common feature is the utilization of the "mass conservation" technique.
As used in cloud resolving or large eddy simulations with a time step of a few seconds, conditionally well-behaved EEBMPC provides a correct numerical solution for the governing differential equations despite the fact that validity of the SM-criterion is never checked.The necessity of its validation is hidden, because a limitation on the time step imposed by dynamics is more restrictive than that imposed by microphysics.Thus, the SM-criterion is always satisfied, and the "mass conservation" technique is never applied.
www.geosci-model-dev.net/5/975/2012/Geosci.Model Dev., 5, 975-987, 2012 Because the SM-criterion determines instantaneous transition between conditionally well-behaved EEBMPC and poorly-behaved EEBMPC, the mechanistic extrapolation of applicability of conditionally well-behaved EEBMPC to regional (global) scales and the utilization for WRF integration time steps on the order of hundreds (thousands) of seconds should be made with caution.As documented in the WRF user guide, limitation on time step permitted for model integration is imposed mainly by dynamics.However, this constraint does not account for additional restrictions due to microphysics imposed by the SM-criterion that ensures explicit Eulerian time integration scheme stability and positivedefiniteness.An occasional violation of the SM-criterion for this time step range determines the necessity of applying the "mass conservation" technique to avoid negativeness of cloud water mixing ratio that makes a conditionally well-behaved EEBMPC become a poorly-behaved EEBMPC (an explicit Eulerian time integration scheme becomes nonstable and non-positive-definite).The eventual feature of a poorly-behaved EEBMPC is that it does not provide a correct numerical solution for the governing differential equations.
Our analysis shows that the source code implementation of single moment (TAO, THOMPSON, and WSM6) schemes and a double-moment MORRISON scheme with a prognostic treatment of precipitating hydrometeors in WRF use the "mass conservation" technique and belong to the conditionally well-behaved EEBMPC class if used for cloudresolving or large-eddy simulations, but they can become poorly-behaved EEBMPC for regional and large scale simulations.
It should be noted that one of the most important aspects of numerical modeling is solving governing differential equations using appropriate numerical methods.If governing differential equations are used, it is obvious that the milestones of applied mathematics, in general, and the theory of finitedifference schemes, in particular, should not be violated.The theory of finite-difference schemes is a branch of science that provides definitions for stability and positive-definiteness of finite-difference schemes (among many others).Our analysis of EEBMPCs in CAM and GFDL AM3 GCM reveals that these extremely important issues are not recognized as essential and crucial.For example, both CAM (Gettelman et al., 2008) and GFDL AM3 GCM (Salzmann et al., 2010) utilize diagnostic equations for precipitating hydrometeors, but the numerical treatment of cloud water remains similar to that used in EEBMPC with prognostic equations, which is discussed above.Additionally, both codes use a mechanistic approach, which is the utilization of equal time substeps because long time steps are used for the host model integration.A feature of these codes is that at a minimum, two substeps are used even if stability and positiveness condition is occasionally satisfied.Moreover, even in a case having a few substeps, stability and positive-definiteness conditions (the SM-criterion) can be violated at different times, altitudes, and spatial locations.
Despite the fact that our analysis is focused on warm rain processes, we highlight that inclusion of the ice phase into consideration makes the SM-criterion even more restrictive because additional solid hydrometeors compete for the available cloud water.For this case the SM-criterion is a necessary but not sufficient condition that should be included in any microphysics scheme whose numerics is based on an explicit Eulerian time integration framework.
As an alternative approach, numerical schemes different from an explicit Eulerian time integration analyzed in this paper can be used in BLK schemes.However, even if the governing differential equations can be solved using stable and positive-definite finite-difference schemes, we would emphasize the need to re-evaluate the validity of the utilization of relatively long time steps in bulk microphysics schemes because of the non-linear dependence of the growth rates of microphysical process on cloud characteristics.It is difficult to expect that the linearization of these growth rates remains valid for periods of time significantly longer than timesteps routinely used in cloud-resolving models.The computational expense of the utilization of smaller timesteps dictated by cloud physics considerations can be prohibitive.However, it is important to develop a numerical framework that is consistent from both the physics and numerics points of view and can be applied to any model regardless of scale.
Our future work is dedicated to the development of a prototype of a so-called open flexible microphysics interface.This interface consists of a suite of different stable positivedefinite time integration schemes (explicit, implicit, and semi-implicit) and contains a process-oriented source code repository with libraries that include functions for calculations of hydrometeor's growth rates due to numerous microphysical process routinely used in different BLK schemes.The distinguishable feature of the openFMI is its ability to "create new bulk microphysics schemes on the fly" eliminating necessity to support the multiple source codes for the BLK schemes in the host model.