Abstract

A benchmark solution that facilitates testing the accuracy, efficiency, and efficacy of moist nonhydrostatic numerical model formulations and assumptions is presented. The solution is created from a special configuration of moist model processes and a specific set of initial conditions. The configuration and initial conditions include: reversible phase changes, no hydrometeor fallout, a neutrally stable base-state environment, and an initial buoyancy perturbation that is identical to the one used to test nonlinearly evolving dry thermals. The results of the moist simulation exhibit many of the properties found in its dry counterpart. Given the similar results, and acceptably small total mass and total energy errors, it is argued that this new moist simulation design can be used as a benchmark to evaluate moist numerical model formulations.

The utility of the benchmark simulation is highlighted by running the case with approximate forms of the governing equations found in the literature. Results of these tests have implications for the formulation of numerical models. For example, it is shown that an equation set that conserves both mass and energy is crucial for obtaining the benchmark solution. Results also suggest that the extra effort required to conserve mass in a numerical model may not lead to significant improvements in results unless energy is also conserved.

1. Introduction

There are many methods for evaluating the performance of a numerical modeling system. For example, the model can be shown to reasonably reproduce an observed meteorological event. Or, it can be demonstrated that the model produces realistic features and evolutions of a phenomenon for which the model was designed to study; for example, a cloud model can be shown to reasonably predict the details of clouds, or a model designed to study the planetary boundary layer can be shown to accurately reproduce the statistical properties of this feature.

The most powerful methods for evaluating a numerical modeling system are comparisons to cases with known results that can be derived through dynamical analysis, or to benchmark solutions that converge under certain conditions. Examples of commonly used analytic and/or benchmark cases include certain mountain wave solutions (e.g., Clark 1977; Dudhia 1993), inertia-gravity waves (Skamarock and Klemp 1994), a nonlinearly evolving cold pool and density current (Straka et al. 1993), and a rising warm thermal (Tripoli 1992; Wicker and Skamarock 1998). Simulations such as these are important for a number of reasons, such as establishing the fidelity of a new numerical modeling system, or testing the accuracy, efficiency, and efficacy of a new numerical technique.

Unfortunately, none of these analytic/benchmark cases include moist processes. Moreover, despite the varying methods used to include moist processes in numerical models, there does not appear to be a commonly agreed-upon method to evaluate a moist model formulation. Typically, a model developer will demonstrate that a model produces reasonable fields of parameters such as vertical velocity, cloud/rainwater mixing ratios, rainfall, etc. for a case of deep, moist convection. While this is an important and necessary step in the development of a numerical model, the lack of a known solution limits the conclusions that can be drawn from such tests. Another common method compares the results from a new model to published results from a different model. While, again, this is an important step in evaluating the fidelity of a new numerical model, this method makes it possible to propagate questionable assumptions through time.

This paper presents a new simulation that can be used as a benchmark for testing numerical models with moisture. The design of the simulation is analogous to the nonlinear warm thermal benchmark case used by Tripoli (1992) and Wicker and Skamarock (1998), but includes phase changes of water vapor and cloud water.

The numerical model used for this study is described in section 2. The dry warm thermal simulation is presented in section 3. Then, the moist test case is presented in section 4. The utility of the test case is demonstrated in section 5 by evaluating common assumptions made in numerical models and the effects these assumptions have on the test case. A summary and conclusions are presented in section 6.

2. The numerical model

a. Governing equations

For the benchmark simulation, the following processes are ignored: hydrometeor fallout, ice-phase microphysics, the Coriolis force, and subgrid-scale turbulence. Under these assumptions, the governing equations1 for a moist atmosphere are

 
formula

(Bannon 2002), where

 
formula

In this paper, Einstein summation convention is used. In (1), the subscript “i = 1, 2, 3” signifies the x, y, and z components, respectively, and δij is the Kronecker delta. The symbols ρa, p, T, rυ, rc, and rt are, respectively, dry air density, pressure, temperature, water vapor mixing ratio, cloud water mixing ratio, and total water mixing ratio. An assumption implicit in the left-hand side of (3) is that the temperature of the liquid water drops is always equal to the air temperature.

In the numerical model, a nondimensional pressure (π) and potential temperature (θ) are solved, which are defined as

 
formula

where p00 = 1000 mb, and

 
formula

To derive governing equations for these two variables, we use the equation of state,

 
p = ρaRT(1 + rυ/ɛ),
(9)

with (2) and (3) to derive a prognostic pressure equation,

 
formula

Then, using (7) and (8), prognostic equations for π and θ are derived, yielding

 
formula

Equations (11) and (12) are used in the present model. It is also worthwhile to point out that a different form of the temperature equation can be derived by using (10) to eliminate the divergence term on the right-hand side of (3), yielding

 
formula

This equation is similar in form to that used in the fifth-generation Pennsylvania State University–National Center for Atmospheric Research Mesoscale Model (MM5; Dudhia 1993).

The thermodynamic equations presented here differ slightly from those used in most numerical models. Traditionally, the specific heats of water vapor and liquid water are ignored in numerical models, so that RmR, cpmlcp, and cvmlcυ, yielding the traditional potential temperature equation,

 
formula

However, it has been shown by several studies (e.g., Lipps and Hemler 1980; Tripoli and Cotton 1981; Bannon 2002) that these terms are necessary for accurate conservation of total energy. Our approach has been to retain all terms in the thermodynamic and pressure equations.

Due to the first term on the right-hand side of (12), the potential temperature equation has a source term even in the absence of phase changes. Thus, potential temperature, as defined by (8), is technically not a conserved variable when water (in any phase) is present. It is interesting to note that an alternative definition for nondimensional pressure given by

 
formula

with a potential temperature defined as

 
formula

results in the following prognostic equations:

 
formula

The final term on the right-hand side of (17) and (18) is nonzero only when phase changes occur,2 since

 
formula

It could be argued that this definition of potential temperature is preferable for use in a numerical modeling system, since it is truly conservative in the absence of moist processes. We have chosen not to use this form, since the formulation of the pressure gradient in the momentum equations becomes more complicated if π̃ is used instead of π; specifically, if the pressure gradient terms were written as ∂π̃/∂xi, additional terms involving the spatial derivatives of rυ and rc would be necessary. We have also found that a model formulated with π̃ and θ̃ is more expensive to run than a model using π and θ, with no perceptible improvements in accuracy.

b. Numerical techniques

The equations actually solved in the numerical model are

 
formula

The advection terms are written as the sum of a flux-form term and a divergence term [i.e., the first two terms on the right-hand side of (20)–(24)]. Writing the advection terms in this way, as opposed to the advective form [second term on the right-hand side of (6)], can sometimes improve conservation of the variable being advected in a numerical model (e.g., Wilhelmson and Chen 1982; Dudhia 1993), although Xue and Lin (2001) point out that it is sometimes possible to formulate an advective form that is numerically equivalent to the flux form with similar conservative properties. However, it should be stressed that not all terms are written in a conservative form, so the numerical model does not exactly conserve mass, momentum, or energy.

In (20), the subscript “0” refers to a hydrostatic base state that varies only in the z direction, and primes refer to deviations from this base state. The hydrostatic equation governing the base state is

 
formula

The density potential temperature, θρ (Emanuel 1994, p. 113), is defined as

 
formula

Following the technique introduced by Klemp and Wilhelmson (1978), the portions of the governing equations that support acoustic waves are integrated on a smaller time step than other terms. The final term on the right-hand side of (20) is a divergence damper, where Kd is a constant and D = ∂uj/∂xj; this term helps maintains stability of the time-splitting technique (Skamarock and Klemp 1992). For the simulations presented here, the model is integrated with third order Runge–Kutta time differencing and fifth-order spatial derivatives for the advection terms (following Wicker and Skamarock 2002).

To account for phase changes, the model uses a saturation adjustment technique, similar to that proposed by Soong and Ogura (1973). In this technique, the equations are advanced forward in two steps: a dynamical step and a microphysical step. In the dynamical step, the model equations are integrated forward with all terms involving phase changes neglected. Then, the microphysics step is applied, in which only the terms involving phase changes are included. This technique is identical to that used by Klemp and Wilhelmson (1978).

Pressure tendencies due to phase changes are included in the microphysics step. Since changes in pressure affect the saturation vapor pressure, an iterative scheme had to be developed for condensation. In this iterative scheme, equations (21)–(24) are advanced forward using a guess for cond. The new values of θ and π are then used to calculate a new value of saturation mixing ratio (rvs), which is then used to calculate a new guess for cond. This cycle is repeated until the newest value of θ converges (to within machine accuracy) to the previous value. During each iteration, the value of cond is determined by the following equation from Rutledge and Hobbs (1983),

 
formula

where Δt is the time step. The iterative technique usually converges in 4–6 iterations.

The ability of the model to conserve mass and energy is evaluated by calculating the total mass (ρt) and the total energy (Et) in the domain during a model integration, where

 
formula

The dry air density, ρa, is determined using the equation of state, (9). The definition of total energy used here is the sum of internal, potential, kinetic, and latent energies—it differs from that used in other studies by the inclusion of the latent energy, which is required to account for the effects of moist processes. For simulations with no flow across the boundaries, and without external sources of mass, momentum, or energy, the domain-total sum of ρt and Et should, theoretically, remain constant over time.

3. The dry simulation

A simulation presented by Wicker and Skamarock (1998) was chosen as the dry reference case. The simulation is two-dimensional, with a domain height of 10 km and width of 20 km. Rigid wall boundary conditions are specified on all four sides of the domain. The initial unperturbed environment is calm (zero initial winds everywhere), hydrostatic, and neutrally stable, defined by a constant potential temperature of 300 K. The initial value of pressure at the surface is 1000 mb, and the initial vertical pressure field is obtained by integrating the hydrostatic equation upwards. A warm perturbation is placed at the center of the domain, which is specified by

 
formula

where

 
formula

xc = 10.0 km, zc = 2.0 km, and xr = zr = 2.0 km. No physical or computational diffusion is applied.

Results of a simulation with 100-m grid spacing after 1000 s of integration are presented in Fig. 1. Similar to the results of Wicker and Skamarock (1998), the thermal rises and expands over time. Two “rotors” develop on the sides of the thermal, while the top of the thermal is stretched. Large θ gradients develop in the middle of the thermal (i.e., within the “arch” spanning between the two rotors).

Fig. 1.

Results of the dry thermal simulation for θ0 = 300 K. (a) Perturbation potential temperature (θ′) is contoured every 0.2 K, and the zero contour is omitted. (b) Vertical velocity is contoured every 2 m s−1, and negative contours are dashed

Fig. 1.

Results of the dry thermal simulation for θ0 = 300 K. (a) Perturbation potential temperature (θ′) is contoured every 0.2 K, and the zero contour is omitted. (b) Vertical velocity is contoured every 2 m s−1, and negative contours are dashed

An interesting extension of this benchmark simulation that, to our knowledge, has not been presented before, is the relative independence of the results to the base-state value of potential temperature. Specifically, if the initial unperturbed state is hydrostatically balanced and neutrally stable, and if an identical initial buoyancy perturbation is applied, the results are similar no matter what value of potential temperature is used. For example, Fig. 2 presents the results of simulations where initial environments of 240 and 270 K are used. In all cases, an identical initial buoyancy perturbation is applied, where buoyancy for the dry case can be defined as

 
formula

For identical initial buoyancy profiles in the two simulations, the specification of initial potential temperature perturbation is given by

 
formula

where θ0 is the reference value of potential temperature for the new simulation, and θ′|300 is the perturbation potential temperature field from the 300-K simulation given by (30). The results displayed in Fig. 2 are very similar to the 300-K reference simulation. Although the fields from these various simulations are not exactly the same, the key point is that the thermal has the same structural details, suggesting a relative independence to the exact thermodynamic value of the neutrally stable base state. Therefore, this property of the solution provides a benchmark for testing new model formulations.

Fig. 2.

Vertical velocity for dry thermal simulations when (a) θ0 = 270 K, and (b) θ0 = 240 K. Contour interval is 2 m s−1; negative contours are dashed

Fig. 2.

Vertical velocity for dry thermal simulations when (a) θ0 = 270 K, and (b) θ0 = 240 K. Contour interval is 2 m s−1; negative contours are dashed

4. The moist simulation

Seeking to obtain a similar result for a moist atmosphere, we again specify the initial environment to be hydrostatic and characterized by exactly neutral stability. In the dry case, it is possible to define neutral stability based on only one thermodynamic variable—potential temperature. However, a moist atmosphere is not as simple. To simplify the specification of the moist base state, two assumptions are made: 1) the total water mixing ratio is constant at all levels, that is, rt = rυ + rc = constant; and 2) phase changes are exactly reversible, that is, Δrυ = Δrvs = −Δrc. Under these two assumptions, a neutrally stable environment can be obtained using one conservative thermodynamic variable (see, e.g., Durran and Klemp 1982). We use the wet equivalent potential temperature, defined for a reversible moist adiabatic atmosphere by

 
formula

(Durran and Klemp 1982; Emanuel 1994, p. 120), where Pd is the partial pressure of dry air. Using (25), (26), and (34), the vertical profiles of π, θ, rυ, and rc can be obtained if values for θe and rt are provided. The value of rt must be greater than rvs at all levels, so that the initial sounding is saturated, and rc > 0 at all levels.

All other parameters are the same as for the dry case, that is, the surface pressure is 1000 mb, the initial wind field is zero, grid spacing is 100 m, and the domain dimensions are as before. No microphysics parameterization is used, other than the assumption of reversible phase changes. Precipitation fallout is not allowed.

The following simple test can be performed using a numerical model to prove that the initial state has been computed accurately, and that the model is configured properly for the benchmark simulation: if no buoyancy perturbation is applied to the model's initial conditions, then no motions should arise. In practice, we have found that small vertical motions (of order 10−4 m s−1) may develop during this test due to truncation errors in the vertical momentum equation (particularly in the balance between the vertical pressure gradient term and the buoyancy term) and/or in the microphysics code (particularly the code that handles condensation), although the present model has been carefully coded so that no vertical motions develop.

For the complete moist thermal simulation, the initial buoyancy field is identical to the dry benchmark case when θ0 is 300 K. For moist conditions, buoyancy is given by

 
formula

Therefore, using (35) for the moist case, (32) for the dry case, and the definition of θρ, the initial θ field is:

 
formula

Since it has been assumed that rt = constant, and rυ = rvs(p, T), the buoyant perturbation has slightly more water vapor and slightly less cloud water than the base-state sounding.

The results after 1000 s of integration for a case in which θe = 320 K and rt = 0.020 are presented in Fig. 3. The results of this moist case are very similar to the results of the dry case (Fig. 1), especially with regards to the structural details such as the two rotors that form on the sides of the thermal and the thin arch that connects them. The moist thermal rises slightly faster than the dry thermal, and after 1000 s the vertical velocity field has higher maximum and minimum values. Nevertheless, the structural details are remarkably similar.

Fig. 3.

Results of the moist thermal simulation for θe = 320 K and rt = 0.020. (a) Perturbation wet equivalent potential temperature (θe) is contoured every 0.5 K. The zero contour is omitted. (b) Vertical velocity, contoured every 2 m s−1

Fig. 3.

Results of the moist thermal simulation for θe = 320 K and rt = 0.020. (a) Perturbation wet equivalent potential temperature (θe) is contoured every 0.5 K. The zero contour is omitted. (b) Vertical velocity, contoured every 2 m s−1

It is important to reiterate that the model formulation for this simulation does not neglect any term in the governing equations. In particular, the specific heat of liquid water (cpl) is included, and the diabatic contribution to the pressure equation is included; it is a common assumption in numerical models to neglect these two effects. Furthermore, the error in total mass and energy conservation is quite small (about 10−4 %), especially compared to model formulations that ignore certain terms in the governing equations (which will be presented in the next section). Given this high degree of accuracy in mass and energy conservation, and the similarity to the dry case, it seems reasonable that this case can be considered a moist benchmark to which moist numerical models can be compared. Additionally, as in the dry case, the simulation proposed here is remarkably insensitive to the values used to define the initial neutrally stable sounding. For example, Fig. 4 shows the results from simulations with θe = 360 K and rt = 0.024, and a case with θe = 280 K and rt = 0.004. Again, the results are nearly identical. These results show that the design of the moist simulation is robust, that is, the correct result is not dependent on a specific initial thermodynamic environment. This is an important point, since it provides further confidence that the results of the simulation truly represent a benchmark solution.

Fig. 4.

Results of moist thermal simulations for: (a) θe = 360 K and rt = 0.024; and (b) θe = 280 K and rt = 0.004. Vertical velocity is contoured every 2 m s−1

Fig. 4.

Results of moist thermal simulations for: (a) θe = 360 K and rt = 0.024; and (b) θe = 280 K and rt = 0.004. Vertical velocity is contoured every 2 m s−1

5. Sensitivity to model formulation

a. Governing equations

The assumptions of reversible phase changes and the absence of hydrometeor fallout clearly make this test case a simplification of reality. However, we have found the case to be valuable for testing the formulation of numerical models, such as defining the governing equations of the model, and for testing numerical techniques involving moist terms. As an example, four different model formulations are tested and presented in this section. Three of these model formulations are found in the literature. The fourth ignores a term that scale analysis suggests has negligible effects on the potential temperature tendency. In all of these cases, only the thermodynamic equation and/or pressure equation is modified.

The four equation sets are summarized in Table 1. Equation set A makes two approximations that are commonly used in nonhydrostatic cloud models: the diabatic contribution to the pressure equation is ignored, and the specific heats of water vapor and liquid water are neglected. This equation set is similar to that used in the Klemp–Wilhelmson model (Klemp and Wilhelmson 1978), in the Penn State–NCAR Mesoscale Model (MM5; Dudhia 1993), and in the Advanced Regional Prediction System (ARPS; Xue et al. 2000), as well as several other numerical models. For equation set B, only the specific heats of water vapor and liquid water are neglected in the thermodynamic and pressure equations. Since the diabatic contribution to the pressure equation is included, this equation set conserves mass. This formulation is similar to that used in the Coupled Ocean–Atmosphere Mesoscale Prediction System (COAMPS; Hodur 1997), and in some respects is similar to models that integrate a density equation rather than a pressure equation (e.g., the Weather Research and Forecasting (WRF) model, Skamarock et al. 2001). In equation set C, the specific heats of water vapor and liquid water are included, as is the diabatic contribution to the pressure equation, but the term involving divergence in the thermodynamic equation [first term on the right-hand side of (12)] is neglected; scale analysis suggests that this term is small and, therefore, negligible. To our knowledge, this equation set has not been used in the literature, but is included here as an example of how approximate forms of the governing equations can be tested numerically. Equation set D uses the ice–liquid water potential temperature (θil) of Tripoli and Cotton (1981) for the thermodynamic equation, where θil is defined as

 
formula

where rl is the liquid water mixing ratio (which is simply rc in the benchmark simulation), and the ice phase has been ignored. Equation set D also neglects the diabatic contribution to the pressure equation, and the specific heats of water vapor and liquid water in the pressure equation. This equation set is similar to that used in Regional Atmospheric Modeling System (RAMS) (Pielke et al. 1992) and in the University of Wisconsin Nonhydrostatic Modeling System (Tripoli 1992).

Table 1.

Summary of thermodynamic and pressure equations

Summary of thermodynamic and pressure equations
Summary of thermodynamic and pressure equations

Before proceeding, it should be noted that all of the equations listed in Table 1 become identical to the “benchmark” equations in the absence of water. Therefore, simulations of the dry thermal case (section 3) are identical for all model formulations presented here. However, if water vapor is present but no liquid water develops (i.e., phase changes either do not occur or are not allowed), then the various equations sets are not identical. This is because the terms involving divergence in the thermodynamic and pressure equations [second terms on the right-hand side of (11) and (12)] only go to zero when water vapor and liquid water mixing ratios are zero.

For these simulations, the initial environment is defined by θe = 320 K and rt = 0.020. The results in Figs. 5 and 6 clearly show the dramatic impact of neglecting terms from the complete thermodynamic and pressure equations—none of the simulations using approximate equations compare well with the benchmark solution (Fig. 3). In all of these cases, the thermal rises much slower than the thermal in the benchmark run, which reaches about 8.2 km. In the θe fields, large undershoots (i.e., anomalously low values, depicted by dashed contours in Fig. 6) develop in all cases.

Fig. 5.

Vertical velocity from moist thermal simulations for various model formulations: (a) equation set A, (b) equation set B, (c) equation set C, and (d) equation set D. See text and Table 1 for details. Contour interval is 2 m s−1

Fig. 5.

Vertical velocity from moist thermal simulations for various model formulations: (a) equation set A, (b) equation set B, (c) equation set C, and (d) equation set D. See text and Table 1 for details. Contour interval is 2 m s−1

Fig. 6.

As in Fig. 5, except for perturbation wet equivalent potential temperature (θe). Contour interval is 0.5 K; negative contours are dashed. The zero contour is omitted. Near the right-hand side of each panel, the height of the top of the thermal is indicated, along with the height of the top of the thermal from the benchmark (BM) simulation

Fig. 6.

As in Fig. 5, except for perturbation wet equivalent potential temperature (θe). Contour interval is 0.5 K; negative contours are dashed. The zero contour is omitted. Near the right-hand side of each panel, the height of the top of the thermal is indicated, along with the height of the top of the thermal from the benchmark (BM) simulation

Other useful conclusions can be drawn from these results. In particular, the w and θe fields from runs A and B are very similar. Both thermals rise to ∼6.9 km, and the vertical motion patterns are nearly identical. This result suggests that the extra effort required to conserve mass in a numerical model (by including the diabatic contribution to the pressure equation) may not lead to significant improvements in results unless total energy is also conserved (as in the benchmark). Despite the similarity of the w and θe fields, the time series of mass and energy errors are dramatically different (Fig. 7). The time series for run A (the short-dashed lines in Fig. 7) has an oscillatory nature, with a period of about 62 s. In contrast, the same time series for run B (dotted lines in Fig. 7) evolves smoothly, with only very minor changes through the simulations. Furthermore, run A loses considerable mass and energy throughout the simulation, with a total mass error that is about 30 times greater than the mass error in the benchmark simulation after 1000 s. Run B has nearly identical mass errors to the benchmark simulation, and has a slight increase in total energy. It is unclear how such dramatically different runs in terms of mass and energy errors can have such similar dynamic and thermodynamic fields (e.g., Figs. 5 and 6).

Fig. 7.

Time series of (a) total mass error, and (b) total energy error from simulations using the five equations sets. Error is expressed as 10−4 % of the total mass/energy at the beginning of the simulation

Fig. 7.

Time series of (a) total mass error, and (b) total energy error from simulations using the five equations sets. Error is expressed as 10−4 % of the total mass/energy at the beginning of the simulation

The results from run C were surprising. Among all the simulations, this run least resembles the benchmark case. The thermal only reached 5.8 km, and the vertical motion pattern is quite different from the other runs. However, it is interesting to note that mass and energy errors from run C are comparable to those in the benchmark simulation (Fig. 7). Apparently, this formulation produces unacceptable results due to an approximation that was made in only one equation, without making a consistent approximation in another equation. Perhaps a “counterbalancing” assumption in the pressure equation would improve the results. It is worthwhile to note that other equations sets make consistent approximations throughout, for example, in equation set B, the specific heats of water are ignored in every equation. On the other hand, equation set A ignores a term in the pressure equation, but does not make a counterbalancing approximation in any other equation, yet the results from equation set A are more acceptable than the results from equation set C. Whatever the reason behind the poor results of equation set C, this test highlights the danger of neglecting terms that may seem unimportant under a scale analysis.

The simulation that used θil as the governing variable (run D) produced w and θe fields that most closely match the benchmark result. The thermal reaches ∼7.6 km, and becomes only slightly distorted in shape. Although the results from the run D formulation are encouraging, this formulation has the largest total mass and total energy errors out of all runs presented here (Fig. 7). On the other hand, θil has other advantages that are not revealed by this test case. As discussed by Tripoli and Cotton (1981), θil is an appropriate variable to use in mixing terms for subgrid-scale turbulence closures. The test case that has been presented here does not include subgrid turbulent processes, nor is a constant background mixing coefficient added. Background mixing was avoided here in order to focus on the formulation of the governing equations, and to highlight the importance of mass and energy conservation in numerical models. It is possible that any disadvantages of θil that are revealed by this simulation are superceded by a more accurate representation of mixing processes.

Although the model used in this paper does not exactly conserve mass, momentum, or energy, the results strongly suggest that conservation of these basic variables can be necessary to obtain accurate results in some instances. This result supports the need to construct numerical models around conservation principles, which is the driving principle behind several recent model development efforts (e.g., Ooyama 2001; Skamarock et al. 2001; Satoh 2002).

The results of this section suggest that the form of the governing equations used in a numerical model may have a profound effect on the simulation of deep, moist convection. On the other hand, one might wonder whether these results only come about due to the unphysical initial environment that must be used to obtain the benchmark solution. Despite the unphysical aspects, this design is required in order for a benchmark solution to be obtained—without this setup, a “correct” solution would not be known, and it would be impossible to objectively evaluate the various model configurations.

Another potential criticism of the simulations presented in this section is the value chosen for rt, which is abnormally high for the imposed temperature sounding. A comparison of simulations with different values for θe and rt reveals that the differences presented here (Figs. 5 and 6) are accentuated over those one would expect to find in more “normal” environments. Nevertheless, it is clear that the mass-conserving and energy-conserving form of the thermodynamic and pressure equations can produce the desired results in all environments, and that these equations should be preferred over approximate equation sets.

We have conducted additional simulations using realistic initial environments (e.g., conditionally unstable and subsaturated initial conditions) to address whether the conclusions drawn from this paper hold for more typical uses of numerical models. We have simulated several forms of deep moist convection, including multicell thunderstorms, supercells, and squall lines using the three-dimensional numerical model with the various equation sets. For these simulations, the Kessler microphysics parameterization is used, hydrometeor fallout is allowed, and a complete subgrid turbulence model is included. Overall, results using the five equation sets are similar. In particular, convective organization is not modified in the any of the tests we have performed (e.g., all of the simulations with environments favoring supercells do produce supercells). However, some of the conclusions noted in this paper are evident in subtle ways. For example, the simulations with the mass- and energy-conserving equation set tends to have the strongest updrafts, the highest cloud tops, and the most rainfall. Simulations with the benchmark equation set tend to have about 10% more rainfall than simulations with equation set A. Based on these results, we have concluded that the form of the governing equations used in a numerical model does have an impact on the results, although perhaps a small impact for most uses.

b. Numerics and assumptions

The moist benchmark simulation proposed here can also be used to test other components of the modeling system. For example, we tested whether saturation adjustment was necessary on the Runge–Kutta time steps. In one case, saturation adjustment was applied once, after the three Runge–Kutta steps were advanced with only dry processes. This simulation was compared to one in which saturation adjustment was applied during each of the three Runge–Kutta steps. The results were nearly identical. Since the former approach is considerably less expensive, we have decided to retain this method in our modeling system.

In another test, we replaced the unapproximated buoyancy term [fourth on the right-hand side of (20)] with the following commonly used approximate formulation for buoyancy:

 
formula

Again, the results were nearly identical to the benchmark simulation, confirming the usefulness of this approximation. These are just a few examples of the utility of the proposed benchmark simulation.

6. Summary and conclusions

A simulation has been designed to evaluate moist nonhydrostatic numerical models. The simulation shares many common characteristics with a dry benchmark simulation that has been used in previous studies. Given a statically stable initial environment, with an identical initial buoyancy perturbation, the moist simulation produces qualitatively similar results to the dry simulation, such as nearly identical structural details. Also, like the dry case, the moist simulation is virtually independent of the exact thermodynamic values of the initial state, provided the environment is neutrally stable. For these reasons, this experimental design is considered to be a benchmark that can be used to evaluate the design of moist nonhydrostatic numerical models.

Using this benchmark, it was shown how the formulation of governing equations can impact the evolution of a rising thermal. Some results obtained from these tests include:

  • Both mass-conserving and energy-conserving equation sets were required to produce acceptable results.

  • When using a common form of the thermodynamic equation, the inclusion of a mass-conserving pressure equation did not improve the result.

  • The neglect of one typically small term from the thermodynamic equation unexpectedly produced the worst results, highlighting the danger of using scale analysis to neglect terms from the governing equations.

  • The results of a simulation using ice–liquid water potential temperature (Tripoli and Cotton 1981) was closest to the benchmark solution, despite having the largest total mass and total energy errors.

The benchmark case is limited for several reasons, including neglect of physical and computational mixing, neglect of ice phase, neglect of hydrometeor fallout, and assumed reversibility of phase changes. Despite these limitations, we have found this case to be extremely useful for evaluating the accuracy and efficiency of numerical techniques and assumptions in a numerical model with moisture.

Acknowledgments

We greatly appreciate the insight provided by Dave Stauffer and Richard James. All figures were created using the Grid Analysis and Display System (GrADS). This work was supported by NSF Grant ATM 9806309.

REFERENCES

REFERENCES
Bannon
,
P. R.
,
2002
:
Theoretical foundations for models of moist convection.
J. Atmos. Sci.
,
59
,
1967
1982
.
Clark
,
T. L.
,
1977
:
A small-scale dynamic model using a terrain-following coordinate transformation.
J. Comput. Phys.
,
24
,
186
215
.
Dudhia
,
J.
,
1993
:
A nonhydrostatic version of the Penn State–NCAR Mesoscale Model: Validation tests and simulation of an Atlantic cyclone and cold front.
Mon. Wea. Rev.
,
121
,
1493
1513
.
Durran
,
D. K.
, and
J. B.
Klemp
,
1982
:
On the effects of moisture on the Brunt–Väisälä frequency.
J. Atmos. Sci.
,
39
,
2152
2158
.
Emanuel
,
K. A.
,
1994
:
Atmospheric Convection.
Oxford University Press, 580 pp
.
Hodur
,
R. M.
,
1997
:
The Naval Research Laboratory's Coupled Ocean/Atmosphere Mesoscale Prediction System (COAMPS).
Mon. Wea. Rev.
,
125
,
1414
1430
.
Klemp
,
J. B.
, and
R. B.
Wilhelmson
,
1978
:
The simulation of three-dimensional convective storm dynamics.
J. Atmos. Sci.
,
35
,
1070
1096
.
Lipps
,
F. B.
, and
R. S.
Hemler
,
1980
:
Another look at the thermodynamic equation for deep convection.
Mon. Wea. Rev.
,
108
,
78
84
.
Ooyama
,
K. V.
,
2001
:
A dynamic and thermodynamic foundation for modeling the moist atmosphere with parameterized microphysics.
J. Atmos. Sci.
,
58
,
2073
2102
.
Pielke
,
R. A.
, and
Coauthors
.
1992
:
A comprehensive meteorological modeling system—RAMS.
Meteor. Atmos. Phys.
,
49
,
69
91
.
Rutledge
,
S. A.
, and
P. V.
Hobbs
,
1983
:
The mesoscale and microscale structure and organization of clouds and precipitation in midlatitude cyclones. VIII: A model for the “seeder-feeder” process in warm-frontal rainbands.
J. Atmos. Sci.
,
40
,
1185
1206
.
Satoh
,
M.
,
2002
:
Conservative scheme for the compressible nonhydrostatic models with the horizontally explicit and vertically implicit time integration scheme.
Mon. Wea. Rev.
,
130
,
1227
1245
.
Skamarock
,
W. C.
, and
J. B.
Klemp
,
1992
:
The stability of time-split numerical methods for the hydrostatic and nonhydrostatic elastic equations.
Mon. Wea. Rev.
,
120
,
2109
2127
.
Skamarock
,
W. C.
, and
J. B.
Klemp
,
1994
:
Efficiency and accuracy of the Klemp–Wilhelmson time-splitting technique.
Mon. Wea. Rev.
,
122
,
2623
2630
.
Skamarock
,
W. C.
,
J. B.
Klemp
, and
J.
Dudhia
,
2001
:
Prototypes for the WRF (Weather Research and Forecasting) model.
Preprints, 14th Conf. on Numerical Weather Prediction, Fort Lauderdale, FL, Amer. Meteor. Soc., J11–J15
.
Soong
,
S-T.
, and
Y.
Ogura
,
1973
:
A comparison between axisymmetric and slab-symmetric cumulus cloud models.
J. Atmos. Sci.
,
30
,
879
893
.
Straka
,
J. M.
,
R. B.
Wilhelmson
,
L. J.
Wicker
,
J. R.
Anderson
, and
K. K.
Droegemeier
,
1993
:
Numerical solutions of a non-linear density current: A benchmark solution and comparisons.
Int. J. Numer. Methods
,
17
,
1
22
.
Tripoli
,
G. J.
,
1992
:
A nonhydrostatic mesoscale model designed to simulate scale interaction.
Mon. Wea. Rev.
,
120
,
1342
1359
.
Tripoli
,
G. J.
, and
W. R.
Cotton
,
1981
:
The use of ice-liquid water potential temperature as a thermodynamic variable in deep atmospheric models.
Mon. Wea. Rev.
,
109
,
1094
1102
.
Wicker
,
L. J.
, and
W. C.
Skamarock
,
1998
:
A time-splitting scheme for the elastic equations incorporating second-order Runge–Kutta time differencing.
Mon. Wea. Rev.
,
126
,
1992
1999
.
Wicker
,
L. J.
, and
W. C.
Skamarock
,
2002
:
Time-splitting methods for elastic models using forward time schemes.
Mon. Wea. Rev.
,
130
,
2088
2097
.
Wilhelmson
,
R. B.
, and
C-S.
Chen
,
1982
:
A simulation of the development of successive cells along a cold outflow boundary.
J. Atmos. Sci.
,
39
,
1466
1483
.
Xue
,
M.
, and
S-J.
Lin
,
2001
:
Numerical equivalence of advection in flux and advective forms and quadratically conservative high-order advection schemes.
Mon. Wea. Rev.
,
129
,
561
565
.
Xue
,
M.
,
K. K.
Droegemeier
, and
V.
Wong
,
2000
:
The Advanced Regional Prediction System (ARPS)—A multiscale nonhydrostatic atmospheric simulation and prediction model. Part I: Model dynamics and verification.
Meteor. Atmos. Phys.
,
75
,
161
193
.

APPENDIX

Definition of Constants and Thermodynamic Variables Not Defined in the Text

  • Variable  Description  Value or definition

  • cp  Specific heat of dry air at constant pressure  1004 J kg−1 K−1

  • cpl  Specific heat of liquid water at constant pressure  4186 J kg−1 K−1

  • cpml  Specific heat of moist air at constant pressure  cpml = cp + cpvrυ + cplrl

  • cpv  Specific heat of water vapor at constant pressure  1885 J kg−1 K−1

  • cυ  Specific heat of dry air at constant volume  717 J kg−1 K−1

  • cvml  Specific heat of moist air at constant volume  cvml = cυ + cvvrυ + cplrl

  • cvv  Specific heat of water vapor at constant volume  1424 J kg−1 K−1

  • g  Acceleration due to gravity  9.81 m s−2

  • Lυ  Latent heat of vaporization  Lυ = Lυ0 − (cplcpv) (TT0)

  • Lυ0  Reference value of Lυ  2.5 × 106 J kg−1

  • rl  Liquid water mixing ratio  rl = rc in this work

  • R  Gas constant of dry air  287 J kg−1 K−1

  • Rm  Gas constant of moist air  Rm = R + Rυrυ

  • Rυ  Gas constant of water vapor  461 J kg−1 K−1

  • T0  Reference temperature  273.15 K

  • ɛ  Ratio of R to Rυ  ɛ = R/Rυ

Footnotes

Corresponding author address: George H. Bryan, 503 Walker Building, University Park, PA 16802. Email: bryan@essc.psu.edu

1

Definitions of terms and notation not presented in the text are provided in the appendix.

2

The tendencies of π̃ and θ̃ will also be nonzero is the presence of hydrometeor fallout, due to the Drl/Dt term.