Abstract

Anelastic and compressible solutions are compared for two moist deep convection benchmarks, a two-dimensional thermal rising in a saturated moist-neutral deep atmosphere, and a three-dimensional supercell formation. In the anelastic model, the pressure applied in the moist thermodynamics comes from either the environmental hydrostatically balanced pressure profile in the standard anelastic model or is combined with nonhydrostatic perturbations from the elliptic pressure solver in the generalized anelastic model. The compressible model applies either an explicit acoustic-mode-resolving scheme requiring short time steps or a novel implicit scheme allowing time steps as large as those used in the anelastic model. The consistency of the unified numerical framework facilitates direct comparisons of results obtained with anelastic and compressible models.

The anelastic and compressible rising thermal solutions agree not only with each other but also with the previously published compressible benchmark solution based on the comprehensive representation of moist dynamics and thermodynamics. In contrast to earlier works focusing on the formulation of moist thermodynamics, the compatibility of the initial conditions is emphasized and its impact on the benchmark solutions is documented. The anelastic and compressible supercell solutions agree well for various versions of anelastic and compressible models even for cloud updrafts reaching 15% of the speed of sound. The nonhydrostatic pressure perturbations turn out to have a negligible impact on the moist dynamics. Numerical and physical details of the simulations, such as the advection scheme, spatial and temporal resolution, or parameters of the subgrid-scale turbulence, have a more significant effect on the solutions than the particular equation system applied.

1. Introduction

Modeling atmospheric dynamics requires numerical methods for solving a system of partial differential equations originating from the Navier–Stokes fluid flow equations. The most complete approach involves the compressible system that is cumbersome because of the presence of fast acoustic modes. These modes are energetically of low relevance for weather and climate, but they impose stringent stability constraints on numerical solvers. Because of that, simplified systems of soundproof equations have been developed in the past (e.g., incompressible Boussinesq, anelastic, or pseudoincompressible) and proved useful in simulations of small-scale and mesoscale low-Mach-number atmospheric flows. Atmospheric regional climate models and general circulation models targeting weather and climate processes at larger scales (synoptic and global) are typically based on an even more simplified system of equations involving hydrostatic balance in the vertical direction (i.e., neglecting fluid accelerations when compared to gravity). Such models have to rely on heavily approximated subgrid-scale parameterizations to represent effects of small-scale processes (turbulence, clouds, precipitation, etc.) on the resolved larger-scale flow. Increasing computational power and the need to better incorporate effects of small-scale processes in weather and climate simulation (clouds and precipitation in particular) provides a strong motivation for the development of models suitable for the entire range of spatial scales, from the microscale to global.

Kurowski et al. (2013, hereinafter KGS13) discussed extension of the soundproof methodologies applied in cloud-scale modeling into larger-scale simulation. It was argued applying a simple scale analysis and documented in numerical simulations that small-scale pressure perturbations (either predicted in the compressible system or diagnosed in the anelastic system) have virtually no effect on the moist thermodynamics. However, computational examples used in KGS13 concerned relatively shallow flows (thermals rising from rest and moist flow over mesoscale topography) and a similar comparison for deep tropospheric flows featuring larger vertical velocities would be desirable. Notwithstanding the scale analysis of Lipps and Hemler (1982), one may wonder if the anelastic framework is suitable for modeling severe convection. The two major issues concern presence of relatively large pressure perturbations that may have a nonnegligible impact on the moist thermodynamics (e.g., saturation adjustment) and an important role of vorticity dynamics in the supercell formation that differ in the anelastic and compressible systems. Moreover, in deep convection, the temperature range encountered by a volume of fluid rising from near the surface to the tropopause is close to 100 K and thus one needs to ensure that the representation of model moist thermodynamics is sufficiently accurate. This issue was emphasized in the idealized test case developed in Bryan and Fritsch (2002, hereinafter BF02).

Keeping in mind our goal to progress toward global modeling, one needs to realize limitations of the fully compressible framework used in KGS13. Grabowski and Smolarkiewicz (2002, their section 5) reported global moist-flow simulations [topographically forced planetary wave and moist Held and Suarez (1994) climate benchmark] applying the anelastic framework with model time steps on the order of 103 s. Since the vertical grid length in those simulations was several hundred meters, applying the explicit compressible model would require model time steps on the order of 1 s to maintain numerical stability with respect to vertically propagating acoustic modes. Smolarkiewicz et al. (2014, hereinafter SKW14) reported a development of a novel implicit (with respect to acoustic, buoyant, and rotational modes) compressible framework that allows application of significantly larger time steps, similar to those used in the anelastic system. An extension of the dry framework of SKW14 to include moist processes is documented in  appendix A and used in simulations of deep convection.

The next section provides a brief overview of the unified fluid-flow modeling system (soundproof/compressible) used in the deep convection simulations and based on the Eulerian–semi-Lagrangian model (EULAG) numerical framework (Prusa et al. 2008; SKW14). We apply two benchmark deep convection simulation cases—namely, BF02 and Weisman and Klemp (1982, hereafter WK82). The former considers two-dimensional deep convection in a fully saturated neutrally stratified environment, and the latter concerns three-dimensional supercell formation in a sheared environment. Model setups and results for the two cases are presented in sections 3 and 4, respectively. Section 5 provides summary and conclusions. In addition,  appendix A summarizes the model numerical framework, and  appendix B discusses idealized simulations that illustrate selected points of the main discussion.

2. The model and model simulations

The research tool employed in this study is the optionally anelastic/compressible EULAG model. The early version of the model was based on the anelastic approximation and proved useful in numerous studies concerning a wide range of low-Mach-number flows (e.g., Prusa et al. 2008). The explicit small-to-mesoscale compressible version of the model originated in Smolarkiewicz and Szmelter (2009) has been recently generalized to moist thermodynamics (KGS13) and extended to an all-scale semi-implicit compressible version facilitating the use of large (e.g., anelastic) time steps (cf. SKW14). A unique feature of the EULAG model is a consistent numerical framework that allows unification of the numerical solver and thus a clear exposition of differences stemming from contrasting mathematical formulations of the governing equations.

The EULAG’s moist thermodynamics is based on a set of prognostic equations for water vapor, cloud condensate, and precipitation mixing ratios (e.g., Grabowski 1998; Grabowski and Smolarkiewicz 2002). In simulations described in section 4, these have been modified to fit warm-rain parameterization applied in WK82. In the anelastic model, the pressure applied in the moist thermodynamics comes from either the environmental hydrostatically balanced pressure profile in the standard anelastic model or is combined with nonhydrostatic perturbations from the elliptic pressure solver as described in KGS13 in the generalized anelastic model. In the compressible model, a prognostic equation for the dry air density is solved and the thermodynamic pressure is subsequently derived from the moist equation of state; see  appendix A for a summary of the unified dynamical core that accounts for moist processes.

Several model configurations are used for the two test simulation setups considered in this paper. These are the standard anelastic model (ANES), the generalized anelastic model with pressure perturbations included in moist thermodynamics (ANEG; see KGS13), the standard anelastic model using compressible time step (ANESc), the implicit compressible model (COMP), the implicit compressible model using the anelastic pressure profile in moist thermodynamics (COMPa), and the explicit compressible model (COMPe). All these are summarized in Table 1.

Table 1.

Description of the EULAG models used in the study with an explanation of main differences in the governing equation sets, model time steps Δt (with a and c subscripts indicating anelastic and acoustic compressible time step sizes, respectively), whether pressure perturbations p′(x, t) are included in or excluded from the moist thermodynamics, and how sound waves are treated in the compressible models.

Description of the EULAG models used in the study with an explanation of main differences in the governing equation sets, model time steps Δt (with a and c subscripts indicating anelastic and acoustic compressible time step sizes, respectively), whether pressure perturbations p′(x, t) are included in or excluded from the moist thermodynamics, and how sound waves are treated in the compressible models.
Description of the EULAG models used in the study with an explanation of main differences in the governing equation sets, model time steps Δt (with a and c subscripts indicating anelastic and acoustic compressible time step sizes, respectively), whether pressure perturbations p′(x, t) are included in or excluded from the moist thermodynamics, and how sound waves are treated in the compressible models.

3. Thermal rising in a cloudy moist-neutral environment

We consider the problem of a moist thermal rising in a cloudy moist-neutral deep atmosphere introduced by BF02. Besides comparing anelastic and compressible solutions, the test provides an opportunity to contrast results obtained with highly accurate representation of moist physics from BF02 (i.e., exact thermodynamic energy and saturated water vapor pressure formulations, comprehensive representation of moisture effects in the governing equations, etc.) with the solutions obtained with simplified thermodynamic energy equation and representation of moist physics as implemented in EULAG [Grabowski and Smolarkiewicz 1996, their Eqs. (1), (3), and (4); Grabowski and Smolarkiewicz 2002].

a. Setup of simulations

BF02 test considers two-dimensional (xz) moist problem limited to reversible adiabatic bulk condensation/evaporation. The domain size is 20 × 10 km2 (x × z) with a uniform grid spacing of 100 m. The atmosphere is neutrally stratified for moist reversible processes; that is, the equivalent potential temperature θe is uniform across the domain. The initial total water mixing ratio has a constant and relatively large value of 20 g kg−1 within the entire domain. With the assumed surface pressure of 1000 hPa, uniform θe dictates that the cloud water of the moist neutrally stratified environment increases from around 7.7 g kg−1 at the bottom boundary to around 19.8 g kg−1 at the model top. Such conditions ensure that the environment remains cloudy when displaced downward by the thermal rise. The atmosphere is initially at rest. The saturated spherical warm thermal (a bubble) is placed 2 km above the ground. The potential temperature perturbation has a 2-K maximum in the center of the bubble and gradually decreases to zero along the 2-km radius. The water vapor mixing ratio is adjusted to maintain water saturation inside the bubble. The solutions are compared after a 1000 s of the thermal rise. The integration time step for the anelastic model is 4 s. The same time step is also used for the implicit compressible model, whereas the explicit compressible (i.e., acoustic) model employs a 0.25-s time step (see Table 1 for a description of the simulations). Further details of the model setup are exactly as in BF02.

Specifics of the uniform-θe condition [i.e., the vertical distribution of the potential temperature θ and water vapor (qυ) and cloud water (qc) mixing ratios that are all used as predicted variables in the EULAG model] depend on particular formulations of the model thermodynamic energy and saturated water vapor mixing ratio. A set of relations between initial/boundary conditions and the governing equations that are necessary and sufficient for their consistency is typically referred to as “compatibility conditions” for partial differential equations (e.g., Temam 2006). In the specific test considered here, these conditions imply that the θ, qυ, and qc profiles corresponding to the uniform θe have to be derived based on the specific model formulation and that applying profiles obtained with a different (i.e., incompatible) formulation may lead to incorrect solutions. To illustrate this very point, we present anelastic and compressible solutions applying θ, qυ, and qc profiles either derived applying the EULAG’s formulation of the moist thermodynamics or as given in BF02—that is, obtained with their comprehensive representation of the moist thermodynamics. Since the conditions at the surface are prescribed, one can compare the temperature and pressure values at the domain top (i.e., at 10 km) obtained with the two formulations. The differences are about 4 K and 2 hPa, and thus one might expect some differences between simulations applying the two sets of environmental profiles as illustrated by model results discussed below.

b. Preamble: Dry anelastic and compressible solutions

In addition to moist simulations, BF02 also presented solutions of a dry thermal rising in a neutrally stratified environment—an analog of their moist case. This allows us to compare dry anelastic and compressible solutions to those shown in BF02.

In general, both anelastic and compressible EULAG solutions proved to be in excellent agreement. The differences between their extreme values at the final time for the vertical velocity w and potential temperature perturbation θ′ were around 1%. The results were also close to the reference solution of BF02 (Fig. 1 therein). The extreme values for w and θ′ fields differed from BF02 by about 4% and 10%, respectively. The latter difference is most likely due to different advection schemes employed in the two models. The EULAG’s nonoscillatory multidimensional positive definite advection transport algorithm (MPDATA) algorithm (e.g., Smolarkiewicz and Margolin 1998; Smolarkiewicz 2006) provides monotone solutions and—in the anelastic system—it is impossible for the minimum (maximum) θ′ to undershoot (overshoot) its initial values of 0 (2) K. The compressible EULAG solution, like the anelastic one, contained neither negative nor larger-than-initial values of θ′, implying that the compressibility of the air did not contribute to the temperature change in both approaches (by design in the anelastic model and because of the low-Mach-number flow regime in the compressible model). Results obtained with the oscillatory MPDATA were closer to those in BF02, with the difference in θ′ extremes reduced from 10% to about 4%. Moreover, the minimum value of θ′ was also negative and not far from the negative value shown in Fig. 1 of BF02. One thus can argue that slightly different θ′ patterns in our dry solutions and in solutions shown in BF02 come from numerical aspects rather than from the effects of air compressibility.

c. Moist solutions

Fields of and w at the BF02’s final time (1000 s) for simulations applying environmental θ, qυ, and qc profiles derived applying EULAG’s formulation of the moist thermodynamics are shown in Figs. 1 and 2, respectively. The figures compare solutions obtained with ANES, ANEG, COMP, and COMPa, all applying the same 4-s time step. In addition, the solutions from the explicit (acoustic) compressible model (COMPe) and the anelastic model (ANESc) both applying a 0.25-s time step are shown. The anelastic ANES and ANEG and compressible COMP and COMPa solutions are almost identical and compare well with Fig. 3 from BF02 in terms of the pattern and vertical reach of the thermal. The consistency of ANES and ANEG (and COMP and COMPa) results clearly indicates that the nonhydrostatic pressure perturbations (~101–102 Pa) are of negligible importance for the latent heating and buoyancy production. This is in agreement with theoretical considerations and numerical simulations discussed in KGS13. The explicit compressible model COMPe provides a solution displaying small-scale features on and w isolines. Results from the ANESc simulation provide an explanation of the main differences between implicit and acoustic compressible moist solutions. Using a smaller time step, the ANESc solution resembles COMPe solution and differs from ANES, ANEG, COMP, and COMPa in the same way as COMPe does. Apparently, reducing time step of the model—and thus increasing the total number of time steps—leads to more frequent modification of the temperature field (owing to more frequent saturation adjustment) and a larger cumulative difference between moist solutions. Similar tendency was also observed in KGS13. It should be stressed that such behavior was not observed in dry simulations from section 3b because dry simulations are less sensitive to the choice of model time step, owing to the absence of complications involved in the advection–condensation problem; see a discussion and computational examples in Grabowski and Smolarkiewicz (1990). In particular, lack of monotonicity for the moist problem even if the monotone advection scheme is used may contribute to the development of small-scale features in the numerical solution.

Fig. 1.

Fields of the equivalent potential temperature perturbations at 1000 s for simulations of thermal rising in a cloudy moist-neutral environment applying various versions of the (left) anelastic and (right) compressible EULAG model (see text for more details). Environmental profiles are prescribed according to the EULAG’s moist physics. Contour interval is shown above the top panels. Model time steps are shown at the bottom of each panel.

Fig. 1.

Fields of the equivalent potential temperature perturbations at 1000 s for simulations of thermal rising in a cloudy moist-neutral environment applying various versions of the (left) anelastic and (right) compressible EULAG model (see text for more details). Environmental profiles are prescribed according to the EULAG’s moist physics. Contour interval is shown above the top panels. Model time steps are shown at the bottom of each panel.

Fig. 2.

As in Fig. 1, but for the vertical velocity.

Fig. 2.

As in Fig. 1, but for the vertical velocity.

Quantitative comparison of model results is shown in Table 2. Because the final shape of the thermal is slightly different between various simulations, the height of the thermal’s center (i.e., the height of the center of mass of the field; Zmean) and the vertical reach of the thermal top (defined as the leading edge height using a 0.1-K threshold; Ztop) are chosen as measures of thermal’s evolution, in addition to the extremes and standard deviations of , w, and p′ and mean values of . The spread of w minima and maxima among various simulations reaches 5%, similarly to the spread of maxima from ANES, ANEG, COMP, and COMPa. The solutions using the compressible time step (i.e., COMPe and ANESc) have about 20% larger maxima. The differences in the solutions seem to be mainly related to the time step rather than to the formulation of the governing equations (soundproof versus compressible). Since the initial minimum is zero, the negative values at the end of simulation have the largest model-to-model variability ranging from −0.16 K for COMP and COMPa to −0.57 K for ANESc and COMPe. In other words, models that use small time step (ANESc and COMPe) have 3–4 times as small minimum values as those using large time steps (ANES, ANEG, COMP, COMPa), with the latter four solutions being the closest to the compressible BF02 solution. The mean values of differ by less than 2% between various EULAG models and are in good agreement with that for BF02. The standard deviation of is about 6% larger for ANES, ANEG, COMP, and COMPa than for BF02 and about 9% larger for the models using the compressible time step. The spread of standard deviation of w between the models does not exceed 6%. The extreme values of the pressure perturbation agree well between BF02 and all EULAG solutions with only 5% difference between the former solution and the compressible models COMP and COMPa. The anelastic models yield slightly more negative values for the minima. For the models using the small time step, those differences are slightly larger.

Table 2.

Comparison of BF02 and anelastic/compressible EULAG solutions for the moist thermal experiment. The columns show the extreme values (min and max) for , w, and p′; the mean values (avg) for ; the standard deviations (std dev) of , w, and p′; the mean height of (Zmean); and the height of the thermal’s leading edge in the middle of the domain (Ztop), all at time of 1000 s.

Comparison of BF02 and anelastic/compressible EULAG solutions for the moist thermal experiment. The columns show the extreme values (min and max) for , w, and p′; the mean values (avg) for ; the standard deviations (std dev) of , w, and p′; the mean height of  (Zmean); and the height of the thermal’s leading edge in the middle of the domain (Ztop), all at time of 1000 s.
Comparison of BF02 and anelastic/compressible EULAG solutions for the moist thermal experiment. The columns show the extreme values (min and max) for , w, and p′; the mean values (avg) for ; the standard deviations (std dev) of , w, and p′; the mean height of  (Zmean); and the height of the thermal’s leading edge in the middle of the domain (Ztop), all at time of 1000 s.

The vertical propagation of the bubble (Zmean and Ztop) is almost the same for all models in Table 2 regardless of the model mathematical formulation or numerical details. This is consistent with the formulation of buoyancy term that is identical in all models, including BF02 [see their Eq. (38)]. The spread between six different EULAG results is at most 1.5% for Zmean (84 m) and Ztop (122 m), with the anelastic solutions achieving slightly lower vertical reach. However, the average difference after traveling the distance of about 4.5 km is smaller than a single 100-m grid length. Overall, it is difficult to find a reliable and systematic difference between anelastic and compressible solutions. Also, modifications in the physics such as including or neglecting pressure perturbations are of a minor importance to the solution (cf. KGS13). Modifications of the model numerics can easily dominate the impact of physical and/or mathematical details of the model formulation. Interestingly, the remarkable agreement between the ANES and COMP solutions is maintained even after 1230 s, when the thermal reaches the domain top and spreads along the upper boundary developing small-scale instabilities (not shown).

Figure 3 presents ANES and COMP model results applying environmental θ, qυ, and qc profiles prescribed as in BF02—that is, assuming BF02’s model formulation to express the constant θe profile through θ, qυ, and qc profiles. The ANES and COMP solutions are almost identical, but they differ significantly from those shown in Figs. 1 and 2 herein and in Fig. 3 in BF02 because the model setup violates the compatibility conditions. In fact, the solutions resemble those from sensitivity simulations discussed in section 5 of BF02 (cf. Figs. 5 and 6 therein). In particular, the solutions show not only positive θe perturbations associated with the imposed initial perturbation but also significant negative perturbations below and above the thermal. These perturbations demonstrate that the environment features nonvanishing θe vertical gradient because only then the displacement of the environmental air (pushed upward above the thermal and dragged below the thermal) can lead to such a pattern. In other words, environmental θ, qυ, and qc profiles prescribed as in BF02 do not result in a uniform θe environment according to EULAG’s moist thermodynamics. Note that this also explains a significantly smaller reach of the thermal (in agreement with results shown in Fig. 6 in BF02) because the environment features nonvanishing stability according to the EULAG’s moist thermodynamics instead of the moist-neutral conditions intended by BF02. As a result, thermal’s rise is suppressed.

Fig. 3.

As in Figs. 1 and 2, but for the initial profiles as given in BF02. Only solutions for ANES and COMP models are shown.

Fig. 3.

As in Figs. 1 and 2, but for the initial profiles as given in BF02. Only solutions for ANES and COMP models are shown.

In summary, results obtained with EULAG’s simplified moist thermodynamics for both anelastic and compressible versions agree between themselves and also with the reference compressible BF02 solution when the environmental θ, qυ, and qc profiles were derived from the uniform θe field applying EULAG’s formulation of the moist thermodynamics. At the same time, the differences between two different compressible numerical models—the one in BF02 and COMP herein—were larger than between anelastic and compressible solutions (ANES and COMP) obtained applying a highly unified numerical framework (see Table 2). When BF02’s environmental θ, qυ, and qc profiles were used, EULAG’s solutions remained close to each other, but were substantially different from BF02 solutions. The latter is explained by the violation of the compatibility conditions in the setup of simulations—that is, using BF02 moist-neutral environmental profiles with the EULAG’s moist thermodynamics. The importance of a consistent derivation of the moist-neutral state for different atmospheric models was also pointed out in Miglietta and Rotunno (2005). Our results are consistent with several other moist deep convection experiments mentioned in BF02 for which no significant differences between various representations of moist thermodynamics were reported. [The role of compatibility conditions (not addressed in BF02) has been clarified by Dr. G. Bryan on his homepage: http://www.mmm.ucar.edu/people/bryan/Code/mbm.html.] A realistic setup of severe convection is discussed in the next section.

4. Idealized simulations of a supercell formation

Kurowski et al. (2011, hereinafter KRZ11) reported a large set of moist deep-convection simulations using observation-based idealized model setups presented in Klemp and Wilhelmson (1978) and WK82. These simulations were conducted using the anelastic version of the EULAG model. One of the conclusions of that study was that simulated storms featured smaller updraft maxima when compared to simulations presented in the original papers. However, such a statement is subject to a few caveats. First, KRZ11 noted a significant sensitivity of the maximum updraft strength to the model spatial resolution (see Figs. 9 and 13 therein), an aspect not addressed in the original studies. Second, the initial sounding in Klemp and Wilhelmson (1978) had only been presented in the graphical form and some uncertainty existed when introducing numerical values for the KRZ11 simulations. Moreover, although analytical formulation of the initial sounding is given in WK82, simulations reported in the previous section suggest that observing the compatibility conditions might require adjustments of the WK82 sounding taking into account EULAG’s formulation of the moist thermodynamics. Finally, there is a possibility that differences in the simulated storm intensity in KRZ11 (as measured by the maximum vertical velocity) come from fundamental differences between anelastic and compressible dynamical frameworks—the issue that KRZ11 were not able to address. We address this issue here by comparing solutions obtained applying various versions of the EULAG model.

To investigate the compatibility conditions, we applied two publicly available codes to calculate the reversible CAPE (i.e., including the condensed water into parcel buoyancy) and the pseudoadiabatic CAPE (i.e., excluding condensed water when calculating the buoyancy) of the WK82 sounding, one written by Dr. George Bryan (http://www.mmm.ucar.edu/people/bryan/Code/getcape.F) and the other by Professor Kerry Emanuel (ftp://texmex.mit.edu/pub/emanuel/BOOK/calcsound.f). The results were compared with a simple adiabatic parcel code based on the EULAG’s thermodynamics. The reversible CAPE was, respectively, 1298, 1242, and 1201 J/kg for the three codes, and the pseudoadiabatic CAPE was, respectively, 1928, 1876, and 2056 J/kg. Assuming that the maximum updraft velocity scales as the square root of CAPE, the differences in the updraft strength implied by the range of CAPE values are on the order of 1 m s−1. This implies that modifying the WK82 sounding to satisfy the compatibility conditions (i.e., to obtain the same CAPE using the EULAG’s formulation of moist thermodynamics) would likely lead to only small changes of the maximum updraft strength. Consequently, adjustments to the sounding as prescribed in WK82 were not considered. The above discussion suggests the BF02 test discussed in the previous section is particularly sensitive to the details of model thermodynamics and that satisfying compatibility conditions in more realistic situations is less significant.

a. Setup of simulations

The simulations follow model setup presented in WK82. Profiles of the relative humidity and potential temperature are defined analytically [cf. Eqs. (1) and (2) in WK82] and represent typical conditions favoring severe deep convection. The setup excludes the Coriolis force and the moist physics is limited to the bulk condensation/evaporation and precipitation processes represented by a Kessler-type warm-rain parameterization following Klemp and Wilhelmson (1978). A unidirectional horizontal velocity profile with a shear layer in the lowest 5 km and the free-tropospheric flow velocity Us = 15 m s−1 was selected (see Fig. 2 in WK82). Anelastic simulations with Us ranging from 0 to 45 m s−1 were discussed in KRZ11. Anelastic and compressible simulations applying different shear strengths were also run in this study, and conclusions were broadly consistent with those for Us = 15 m s−1. To keep the solution in the center of the domain, a Galilean transformation with a constant horizontal velocity of 9.15 m s−1 was applied. The domain size (length × width × height) was 128 × 128 × 17.5 km3 with horizontal and vertical grid lengths of 2 km and 350 m, respectively. The integration time was 120 min with time steps of 10 s for anelastic and implicit compressible simulations and 0.5 s for the explicit compressible model (see Table 1 for a description of the models).

The horizontal domain size applied in current simulations is significantly larger than the one in WK82. The main reason comes from different lateral boundary conditions: instead of the open boundaries as in WK82, current simulations apply periodic conditions and a sponge zone near lateral boundaries. The sponge zone (Davies 1976) extends over six model grid points and features relaxation toward the prescribed profiles for all model variables with the time-scale of 50 s. The main difference between open and periodic lateral boundaries is that the horizontally averaged vertical velocity at all levels has to vanish when periodic lateral boundary conditions are used.1 With open lateral boundaries, the model is capable in developing the mean ascent across the domain associated with convectively-driven low-level convergence and upper-level divergence. KRZ11 discussed the impact of lateral boundary conditions in their simulations by comparing results obtained with increasingly large horizontal computational domains, up to 256 × 256 km2. The domain applied here is sufficiently large so the impact of the lateral boundary conditions can be neglected. Rigid-lid free-slip boundary conditions are applied at lower and upper boundaries, as in WK82 and KRZ11. As in the previous section, six model configurations are considered (Table 1), namely ANES, ANEG, COMP, COMPa, COMPe and ANESc.

b. Results

Overall, the results applying different model configurations are in excellent agreement as illustrated by the following discussion. Figures 4 and 5 combine various model fields to highlight the agreement. Before discussing the figures, a comment concerning the pressure field is in order. The pressure perturbation fields shown in the figures exclude the horizontally averaged quasi-hydrostatic pressure profile that develops in a response to the mean (horizontally averaged) buoyancy perturbation. In the compressible model, the horizontally averaged pressure profile also includes a constant component that forms because of the domain-averaged temperature increase. We will refer to this time-evolving component as the compressible pressure offset. We include an extensive discussion of various aspects of the anelastic and compressible pressure fields in  appendix B. It provides computational examples illustrating development of the pressure offset in the compressible model and formation of hydrostatically balanced pressure profiles in the anelastic and compressible models. We stress that these components are removed from the pressure field in the figures (as only nonhydrostatic pressure perturbations contribute to the vertical momentum equation) to document similarities between model solutions. However, they are included in all model computations, including the moist thermodynamics (unless purposely omitted). Nonetheless, their magnitude is so small (101–102 Pa) that their effect on the moist processes is insignificant, as illustrated by the discussion below.

Fig. 4.

Comparison of four different EULAG solutions (ANES, ANEG, COMP, and COMPa) after 120 min from the WK82 storm-splitting experiment for the Us = 15 m s−1 case. Each panel represents a horizontal cross section through the supercell and shows the surface rain mixing ratio (light blue lines; contours at 0.1, 1.5, 3, 4.5 g kg−1), surface cold pool edge (thick dark blue line; defined by the θ′ = −0.5-K contour), updraft strength at 4900 m (black solid lines; contour interval of 4 m s−1), and pressure perturbations at 2800 m (red lines; contour interval of 15 Pa). Positive (negative) values are marked with solid (dashed) lines. Arrows show surface horizontal flow. As in WK82, the mean flow of 12 m s−1 was subtracted from the x component of horizontal velocity. Only half of the solution is shown for each model. The horizontal grid length is 2000 m.

Fig. 4.

Comparison of four different EULAG solutions (ANES, ANEG, COMP, and COMPa) after 120 min from the WK82 storm-splitting experiment for the Us = 15 m s−1 case. Each panel represents a horizontal cross section through the supercell and shows the surface rain mixing ratio (light blue lines; contours at 0.1, 1.5, 3, 4.5 g kg−1), surface cold pool edge (thick dark blue line; defined by the θ′ = −0.5-K contour), updraft strength at 4900 m (black solid lines; contour interval of 4 m s−1), and pressure perturbations at 2800 m (red lines; contour interval of 15 Pa). Positive (negative) values are marked with solid (dashed) lines. Arrows show surface horizontal flow. As in WK82, the mean flow of 12 m s−1 was subtracted from the x component of horizontal velocity. Only half of the solution is shown for each model. The horizontal grid length is 2000 m.

Fig. 5.

As in Fig. 4, but for the vertical cross sections (yz) through the convective cores for ANES and COMP solutions at (a) 60 and (b) 120 min. Black contours are vertical velocity every 5 m s−1, shaded regions show presence of rainwater with qr > 1 g kg−1, light blue isolines mark the cloud edge (qc > 0.1 g kg−1), dark blue thick line indicates location of the near-ground cold pool (θ′ = −0.5 K), and red lines are pressure perturbations contoured with a 40-Pa interval; negative values are dashed. Only half of the solution is shown for each model. The vertical grid length is 350 m.

Fig. 5.

As in Fig. 4, but for the vertical cross sections (yz) through the convective cores for ANES and COMP solutions at (a) 60 and (b) 120 min. Black contours are vertical velocity every 5 m s−1, shaded regions show presence of rainwater with qr > 1 g kg−1, light blue isolines mark the cloud edge (qc > 0.1 g kg−1), dark blue thick line indicates location of the near-ground cold pool (θ′ = −0.5 K), and red lines are pressure perturbations contoured with a 40-Pa interval; negative values are dashed. Only half of the solution is shown for each model. The vertical grid length is 350 m.

Figure 4 compares anelastic and compressible supercell solutions at 120 min in a manner similar to WK82. The plots represent horizontal cross sections through the simulated structure at 120 min and document the storm-relative horizontal flow near the ground, the extent of the surface cold pool, the surface precipitation rate, the vertical velocity at 4900 m, and the nonhydrostatic pressure perturbations at 2800 m. Only half of the domain is shown (i.e., the left cell for ANES and COMP and the right cell for ANEG and COMPa) because the solutions remain symmetric with respect to the y = 0 axis. As the figure documents, all model configurations reproduce storm splitting in a consistent way, and the solutions differ only in minor details. The location of convective cores is similar in all cases, so are the low-level winds and the updraft intensities. Rain formation and fallout progress similarly and the surface precipitation areas are almost identical in all solutions. The shape of the cold pool formed as the negatively buoyant precipitation-laden air hits the surface and spreads horizontally seems to be almost unaffected by the model configuration. The anelastic (diagnosed) and compressible (prognosed) pressure perturbations compare well in terms of the spatial distribution and the magnitude. The largest differences (about 30 Pa) are ahead of the storm (i.e., to the right of the cold pool edge) and between the two convective cells (i.e., near the y = 0 axis).

Figure 5 shows yz cross sections of ANES and COMP solutions at 60 and 120 min, at the same x location for both simulations, close to the maximum vertical velocity. As in Fig. 4, we take advantage of the solution symmetry and only show either right or left half of the solution. Similar to Fig. 4, the figure documents strong similarities between the two solutions. The pattern and magnitude of all fields, including the nonhydrostatic pressure perturbations, also compare well between the models.

Small differences between solutions obtained with the six model configurations are further supported by the evolutions of the vertical velocity and cloud and rainwater mixing ratio maxima shown in Fig. 6. The w maxima from different simulations are close for all six model configurations with the solution envelope in the first hour defined by ANESc (the lowest values) and COMP (the highest values). In general, compressible models yield a few-percent-larger vertical velocity maximum wmax than do anelastic models for the same time step. This difference is attributed to the baroclinic production of vorticity that is abbreviated in the anelastic model.2 As it will be shown in the sensitivity experiments, the difference slightly increases for larger updraft velocities. This is the first systematic difference between the two mathematical systems observed in the study. Within each group (i.e., anelastic versus compressible), sensitivity to the time step is also similar, with the largest maxima for the largest time steps. All realizations of convection for a given time step are similar to each other, regardless of the model configuration, but are slightly different for different time steps. Insignificant differences between ANES and ANEG as well as COMP and COMPa solutions imply that not only compressible and anelastic dynamics differ little in this case, but also that including or excluding pressure perturbations into the moist thermodynamics has virtually no impact. These differences are much smaller than those resulting from various modifications of the model numerics (e.g., the time step as shown in Fig. 6; advection scheme), domain geometry (e.g., size and horizontal/vertical grid length as documented in KRZ11), or subgrid-scale parameterizations (as will be illustrated in section 4c). The evolutions of the cloud water and rainwater mixing ratio maxima are practically indistinguishable between models applying either large or small time steps. Extending simulations up to 240 min (i.e., beyond those documented in WK82 and KRZ11) continues to show excellent agreement between the simulations.

Fig. 6.

Evolutions of the (top) maximum vertical velocity, (middle) maximum cloud water mixing ratio, and (bottom) maximum rainwater mixing ratio for ANES, ANEG, ANESc, COMP, COMPa, and COMPe supercell simulations. The insert in the top panel shows detailed evolution of wmax for times between 20 and 40 min.

Fig. 6.

Evolutions of the (top) maximum vertical velocity, (middle) maximum cloud water mixing ratio, and (bottom) maximum rainwater mixing ratio for ANES, ANEG, ANESc, COMP, COMPa, and COMPe supercell simulations. The insert in the top panel shows detailed evolution of wmax for times between 20 and 40 min.

Small differences in the pressure field shown in Figs. 4 and 5 are apparently insignificant for the model dynamics as documented in Fig. 6. This is most likely because not the pressure magnitude but the pressure gradient is relevant for the dynamics (and the pressure gradients compare well; cf. Figs. 4 and 5), or alternatively because the magnitude of the differences (typically a fraction of 1 hPa) is simply too small to have any appreciable effect. Moreover, the differences are also insignificant for the model thermodynamics despite the fact that in this case the pressure magnitude does matter. Nevertheless, the differences deserve attention to better understand physical processes in anelastic and compressible models. Figure 7 shows profiles of horizontally averaged pressure perturbation fields from the initial hydrostatically balanced profile. For the ANES model, this is the horizontally averaged output from the elliptic pressure solver. For the compressible model, the pressure perturbation also excludes the time-evolving volume-averaged offset that develops as the simulation progresses. Evolution of this volume-averaged mean pressure perturbation is shown in Fig. 7c. This systematic offset develops in the compressible model because of the continuous latent heat release and gradual increase of the mean temperature of the domain. Since we apply periodic lateral boundary conditions and rigid-lid lower and upper boundary conditions, the simulations represent processes inside a closed box, and the increase of the mean temperature has to lead to the increase of the mean pressure. This effect is absent in the anelastic system because the total pressure combines the initial pressure profile and dynamics-driven pressure perturbations from the elliptic solver. The development of the pressure offset in the compressible system and the hydrostatic pressure profile in the anelastic and compressible systems are discussed in detail in  appendix B.

Fig. 7.

Profiles of horizontally averaged pressure perturbations for (a) ANES and (b) COMP models. (c) Evolution of the pressure offset for the COMP model that was subtracted from profiles shown in (b).

Fig. 7.

Profiles of horizontally averaged pressure perturbations for (a) ANES and (b) COMP models. (c) Evolution of the pressure offset for the COMP model that was subtracted from profiles shown in (b).

c. Sensitivity simulations

Sensitivity simulations in KRZ11 showed a considerable dependency of wmax for the entire simulation on the horizontal and vertical resolution. The greatest value of wmax increased from around 29 to 39 m s−1 when the model grid lengths changed from Δxz = 2000/350 m to Δxz = 1000/175 m. We conducted additional simulations using anelastic and compressible models with grid lengths of Δxz = 1000/175 m. The results proved similar to KRZ11, with large sensitivity to model grid length for both compressible and anelastic models. Perhaps more importantly, with much stronger updrafts for the finer resolution, the difference between cloud water and vertical velocity maxima between anelastic and compressible solutions remained as small as in Fig. 6.

To further document similarities and differences between anelastic and compressible solutions, we also considered the impact of the mixing length Λ in the EULAG’s subgrid-scale scheme based on the prognostic equation for the turbulent kinetic energy (Margolin et al. 1999). The mixing length is typically selected considering the model grid length. This is straightforward when each spatial direction has the same grid spacing. However, when horizontal and vertical grid lengths differ significantly, selecting Λ is more ambiguous. Figure 8 shows wmax for different Λ values for ANES and COMP solutions. Since model grid lengths differ significantly (Δxz ≈ 6), the range of Λ is used to demonstrate the sensitivity. The maximum updraft velocity increases from about 26 to 47 m s−1 for Λ decreasing from Δx to zero. The anelastic and compressible models provide consistent solutions with only a few-percent difference for each Λ. The difference increases for larger updrafts, which is the first systematic disparity between the two approaches. Nonetheless, even for the largest updraft velocities reaching almost 50 m s−1 (Mach number 0.15) the anelastic and compressible solutions remain close to each other, with only about 4%-weaker updrafts for the anelastic model. The impact of Λ on the supercell’s spatial structure is shown in Fig. 9. Apparently, the subgrid-scale transport controls the formation of convective structures and affects formation of the cold pool. The differences for different Λ values are significantly larger than those resulting from different mathematical formulations (cf. Fig. 4).

Fig. 8.

Maximum vertical velocity for the entire simulation for ANES and COMP simulations as a function of the mixing length in the subgrid-scale turbulence scheme. Implicit large-eddy simulation represents simulation without explicit subgrid-scale model.

Fig. 8.

Maximum vertical velocity for the entire simulation for ANES and COMP simulations as a function of the mixing length in the subgrid-scale turbulence scheme. Implicit large-eddy simulation represents simulation without explicit subgrid-scale model.

Fig. 9.

As in Fig. 4, but for ANES model applying four different values of the mixing length.

Fig. 9.

As in Fig. 4, but for ANES model applying four different values of the mixing length.

5. Summary and conclusions

After exploring shallow small-scale and mesoscale anelastic and compressible flows in KGS13, the focus of the current study is on the simulation of deep convection. The temperature range encountered by a parcel rising from near the surface to the tropopause and the magnitude of vertical velocities in deep convection justifies a closer look at numerical solutions obtained applying the two contrasting frameworks. The overall goal of this research is to develop universal and efficient modeling approaches valid across the entire range of scales, from the small-scale dynamics to planetary flows. As argued in KGS13, compressible dynamics provides the most comprehensive framework for the all-scale atmospheric modeling but impose computational limitations that are difficult to overcome. Nevertheless, compressible dynamics can provide reference solutions for simplified frameworks, such as anelastic dynamics. In the case of moist processes, additional issue is the incorporation of pressure perturbations into the moist thermodynamics (cf. KGS13).

We presented a detailed comparison between anelastic and compressible solutions for the two deep convection benchmark cases: a two-dimensional thermal rising in a deep moist-neutral saturated environment (BF02) and a three-dimensional severe convection leading to the supercell formation (WK82). Numerical experiments were performed using consistent numerical framework of the EULAG model, which enables a better control of numerical effects and consequently a more credible exposition of the differences resulting from contrasting formulations of the governing equations. The implicit compressible model permits the use of large time steps (i.e., the same as in the anelastic model in simulations presented here) and mitigates the severe limitation of the compressible system.

Our simulations did not confirm the high sensitivity of model results to the details of the moist thermodynamics formulation as suggested in BF02. In contrast, the anelastic and compressible EULAG solutions based on EULAG’s simplified thermodynamics agreed not only with each other, but also with the compressible benchmark solution from BF02—the latter based on the comprehensive representation of moist thermodynamics. The explanation of sensitivity reported in BF02 is in violation of the compatibility conditions, specific for each set of the governing equations, which is an aspect not addressed in BF02. A base state that is moist neutral under one formulation of moist thermodynamics may have nonzero moist stability when the governing equations change—for example, through the neglect of water vapor or liquid water contributions to the specific heat.

As far as WK82 benchmark is concerned, the anelastic and compressible models provided consistent solutions as well. The patterns of the vertical velocity, potential temperature, and moisture variables were almost identical in all simulations. The evolutions of the maximum vertical velocity were close to each other, with slightly stronger updrafts for the compressible system. This difference, however, was insignificant for the evolution of the supercell. The pressure perturbations, diagnosed/prognosed in the anelastic/compressible framework, also compared well, except for the buildup of the quasi-hydrostatic component due to the latent heat release in a periodic computational domain. The resulting mean temperature increase leads to the development of a pressure offset in the compressible model. In contrast, the elliptic pressure solver of the anelastic model derives the pressure perturbations without the development of any mean component (see the appendix in KGS13). This difference results from the choice of boundary conditions and can be argued to diminish for open conditions that allow propagation of acoustic modes out of the domain. The quasi-hydrostatic component along with the nonhydrostatic pressure perturbations were found to have a negligible impact on the saturation adjustment even for the updrafts exceeding 30 m s−1, consistent with findings of Wilhelmson and Ogura (1972). Perhaps most importantly, differences in the numerical/physical details (such as the advection scheme or details of the subgrid-scale transport) turned out to have a larger impact than those resulting from application of either of the two frameworks.

Comparison of numerical efficiency of the model versions applied in this study shows that the anelastic and implicit compressible codes, both using a large anelastic time step dta, require about the same computing time and are roughly 0.5dta/dtc faster than the explicit compressible code, where dtc is the acoustic compressible time step. Therefore, a benefit of using larger time steps in anelastic and implicit compressible models is partially reduced by the algorithms employed to solve the Poisson or Helmholtz equations, respectively. As argued in the introduction, the advantage of using large time steps becomes increasingly important for larger-scale flows. This paves the way for computationally efficient numerical simulation at synoptic and global scales, where vorticity dynamics differ significantly between anelastic and compressible frameworks. This issue is currently being addressed and will be reported in forthcoming publications.

Acknowledgments

This material is based upon work partially supported by the U.S. Department of Energy under Award DE-SC0006748. The views and opinions of the authors expressed herein do not necessarily state or reflect those of the U.S. Government or any agency thereof. WWG was also partially supported by the NSF Science and Technology Center for Multiscale Modeling of Atmospheric Processes (CMMAP), managed by Colorado State University under Cooperative Agreement ATM-0425247. PKS acknowledges support by funding received from the European Research Council under the European Union’s Seventh Framework Programme (FP7/2012/ERC Grant Agreement 320375). Dr. Zbigniew Piotrowski advanced massively parallel code of the model used in simulations described in this paper under partial support of the project “Towards peta-scale numerical weather prediction for Europe” implemented within the Homing Plus programme of the Foundation for Polish Science and cofinanced by the Regional Development Fund of the European Union.

APPENDIX A

A Unified Anelastic/Compressible Framework for Moist Flows

A numerical framework for consistent integrations of soundproof and the fully compressible nonhydrostatic equations of motion for inviscid, adiabatic, all-scale dry atmospheric flows has been documented in SKW14. Because the design of the compressible solvers derives from the preexisting anelastic solvers, viscous, diabatic, and moist processes can be easily incorporated into the compressible integrators; in particular, without affecting the moist anelastic schemes documented in Grabowski and Smolarkiewicz (1990, 1996, 2002). However, these processes enter the dynamical core, and the elliptic Helmholtz solver in particular, as additional forcings. Here we highlight this aspect of the framework, limiting the mathematical formalism of SKW14 to that relevant for the current paper. Thus, we restrict the presentation to the anelastic and compressible PDEs, and write the governing equations in a physically intuitive nonrotating Cartesian vector form, in abstraction from the model geometry and the coordinate frame adopted. Furthermore, in the discussion that follows we do not include equations representing evolution of the moist phase but only account for their contributions to the dynamical core.

The generalized anelastic/compressible dynamical core that accounts for moist, diabatic, and viscous forcings can be compactly written as

 
formula
 
formula
 
formula

Here, vector u denotes the flow velocity, θ is the potential temperature, while φ and ϱ denote a perturbation of the generalized Exner pressure and generalized density of dry air. Furthermore, the differential operators of the total derivative d/dt and the nabla = (∂x, ∂y, ∂z) have their generic meaning; so, d/dt = ∂/∂t + u · . On the rhs of the momentum equation [Eq. (A1)], FE combines all explicit forcings, including moist part of the buoyancy (cf. Grabowski and Smolarkiewicz 2002) and the divergence of the viscous stress.A1 On the rhs of the entropy equation [Eq. (A2)], the source–sink term H subsumes all diabatic effects including latent heating and diffusion. The subscript b indicates the basic (reference) hydrostatically balanced state of the anelastic asymptotic expansion, whereas primes denote perturbations with respect to a balanced ambient state θa, πa—compatible with the governing PDEs—that can, but does not need to, coincide with the basic state θb, πb. Standardly, g = (0, 0, −g) denotes the gravitational acceleration. The generalized density and pressure variables ϱ and φ are defined, respectively, for the [anelastic, compressible] PDE sets as

 
formula

together with corresponding dimensionless coefficients

 
formula

where cp is the specific heat at constant pressure, θ0 is a constant reference value of θ, qυ is water vapor mixing ratio, qt denotes the total water mixing ratio, and ϵ = Rd/Rυ. While the generalized system [Eqs. (A1)(A3)] takes a form of compressible equations, its interpretation depends on the definition of the generalized density used in Eq. (A3), either as a prescribed problem parameter for the anelastic system [Eq. (A4)] or as a dependent prognostic variable with the associated constitutive law

 
formula

where ξ:= Rd/(cpRd) and p0 is a constant reference pressure.

Combining ϱ·[Eq. (A1)] with u·[Eq. (A3)], and ϱ·[Eq. (A2)] with θ′·[Eq. (A3)], and combining the rhs of Eq. (A3) with the total derivative dϱ/dt on the lhs, leads to the set of conservation laws

 
formula
 
formula
 
formula

where Ru and Rθ symbolize the rhs of Eqs. (A1) and (A2), respectively, and ⊗ denotes the tensor product. The equations for momentum [Eq. (A7)] and potential temperature [Eq. (A8)] can be viewed as the generic conservation law

 
formula

consistent with the Lagrangian form

 
formula

for Eqs. (A7) and (A8), respectively; here, ψ—a specific variable expressed per unit of mass—symbolizes the three components of the velocity vector (viz. specific momenta) and potential temperature perturbation (tantamount of specific dry entropy), while denotes the associated rhs. Importantly, the generalized dry air mass continuity equation [Eq. (A9)] is a special case of Eq. (A10), with predetermined ψ ≡ 1 and ≡ 0 for all (x, t) and both sets of governing PDEs. This makes the mass continuity distinct from conservation laws for specific dependent variables ψ, and its role for the design of the framework is summarized next.

A standard template for the EULAG’s second-order anelastic integrators takes compact symbolic form

 
formula

where i, tn refers to the position on a collocated grid discretizing an (x, t) model domain. The transport operator symbolizes a second-order-accurate nonoscillatory advection scheme MPDATA (cf. SKW14 and references therein) defined in Eq. (A12) in terms of its entries: the auxiliary field , an O(δt2)-accurate estimate of the transportive momentum Vn+1/2 = ϱun+1/2 at an intermediate time level tn+1/2 = tn + 0.5δt, and the prescribed anelastic density ϱ = ρb(z) in Eq. (A4). A related template for compressible equations takes the form

 
formula

assuming the preceding integration of the prognostic mass continuity equation in Eq. (A9)

 
formula

where (x, t) ≡ 1. The integration in Eq. (A14) supplies Eq. (A13) with both the updated value of ϱn+1 and the Vn+1/2 defined as the cumulative mass flux in the solution [Eq. (A14)]. The latter assures compatibility of the transport operator in Eq. (A13) with the mass continuity; for example, it preserves a local constancy of ψ in conservative advection, in analogy to the anelastic integrator [Eq. (A12)].

Given the template (A13), an iterative acoustic algorithm accounting for the nonlinearity of the pressure gradient force on the rhs of the momentum equations can be compactly written as

 
formula

where and are shorthand for applied to and in (A13) and subsumed with contributions from the explicit forcings 0.5δtHn+1 and , respectively. Furthermore,

 
formula
 
formula

where the availability of updated water mixing ratios, such as in Eq. (A16), is ensured by design of the moist scheme (Grabowski and Smolarkiewicz 2002). Throughout Eqs. (A15)(A17), the index ν = 1, …, Nν numbers the iterations, with the first guess generated by advecting full θ,

 
formula

With this design, the solution is fully second-order accurate even for Nν = 1, and Nν = 2 gives already close approximation to the trapezoidal integral (Smolarkiewicz and Szmelter 2009); relevant calculations reported in this paper use Nν = 3.

The scheme outlined in Eqs. (A15)(A17) contains a fully implicit trapezoidal integral of the dry buoyancy, whereas pressure perturbations (viz. acoustic modes) and the coefficient Θ depending on full potential temperature and water mixing ratios are integrated explicitly. To derive the closed-form expression for the velocity update, we substitute the potential temperature perturbation in the buoyancy term of the momentum equation with the rhs of the entropy scheme and gather all terms depending on uν on the lhs of the momentum scheme while dropping the spatial grid index i everywhere, as all dependent variables, coefficients, and terms are collocated in Eqs. (A15)(A17). This results in a system of three linear algebraic equations with three unknown components of the velocity vector uν at each point of the collocated grid:

 
formula

the closed-form solution of which may be compactly written as

 
formula

where and = L−10.5δtΘν−1 denote a 3 × 3 matrix of known coefficients (cf. SKW14 and references therein). In each iteration ν the velocity update in Eq. (A20) uses the thermodynamic pressure in Eq. (A16), and the total potential temperature (required in the coefficient Θ) gets updated according to Eq. (A17). The potential temperature perturbation θ′ is updated according to Eq. (A15), upon completion of the velocity update for ν = Nν.

The acoustic scheme summarized above is an adaptation of the soundproof pseudoincompressible algorithm (Smolarkiewicz and Dörnbrack 2008)—an extension of the anelastic algorithm (Grabowski and Smolarkiewicz 2002; Prusa et al. 2008). In essence, the acoustic scheme replaces the anelastic density and the elliptic pressure perturbation with the prognosed thermodynamic density and pressure. Conversely, with the density prescribed as ϱ = ρb(z), and thus the mass continuity [Eq. (A9)] reduced to

 
formula

the anelastic algorithm consistent with the acoustic scheme solves the Poisson equation implied by Eq. (A21) applied to Eq. (A20):

 
formula

This consistency minimizes numerical differences between the anelastic and compressible algorithms and provides the reference solution for large-time-step compressible schemes for atmospheric flows. Furthermore, the acoustic scheme forms the foundation of the semi-implicit compressible models with large (soundproof) time steps, discussed next.

To derive the large-time-step semi-implicit compressible integrators, SKW14 employ the evolutionary form of the gas law [Eq. (A6)], rather than the gas law itself [Eq. (A16)]. Combining ϱd/dt[Eq. (A6)] with φ·[Eq. (A3)] leads to the conservation-law form [Eq. (A10)] for φ:

 
formula

with the rhs forcing

 
formula

where ϕa = cpθ0πa(x), ϕ = φ + ϕa, ξ was defined in Eq. (A6), and the explicit forcing Π is given as

 
formula

with Qυ symbolizing the rhs forcing of the qυ equation. The resulting equation [Eq. (A23)] is integrated to O(δt2) using the template

 
formula

with the auxiliary field under the transport operator defined as and the implicit forcing . Integrating the pressure equation to the first-order accuracy in time still suffices for the second-order-accurate model solution, as φ enters Eq. (A20) with the δt factor. Denoting the first term on the rhs of Eq. (A26) shortly as and regrouping all terms on the rhs leads to the implicit Helmholtz problem for φν

 
formula

to replace the explicit thermodynamic pressure in Eq. (A16) of the acoustic scheme. The resulting Helmholtz operator is composed of the three Poisson operators of the form in Eq. (A22), each using a different effective density—1, ϱn+1ϕa, and ϱn+1, respectively—and entering the problem with different weights. The second and the third Poisson operators combine into the term proportional to u · ϕa, so the expression in Eq. (A27) could be simplified. However, the form of Eq. (A27) is beneficial, as it allows us to build the Helmholtz solver from the Poisson solver with minimal changes to the EULAG model code.

APPENDIX B

Effects of Heat Sources on the Mean Pressure in Anelastic and Compressible Models

To investigate the development of the mean pressure offset in the compressible system, and to highlight the role of lateral boundary conditions in this problem, we developed a simple one-dimensional compressible model that mimics physical processes taking place in multidimensional compressible simulations. Details of the model and results of simulations illustrating the development of the pressure offset are presented in the next section. The main conclusion is that the offset develops in simulations with periodic lateral boundary conditions (i.e., for a closed box) and that the mean pressure should not change if more natural truly open lateral boundary conditions are used. Subsequently, we consider an idealized problem of horizontally uniform constant heating applied within a layer in the lower troposphere in the two-dimensional nonrotating stratified dry atmosphere. The simulations document hydrostatic adjustment in the anelastic and compressible models in response to the continuous heating and aid interpretation of supercell simulations discussed in section 4.

a. Development of the pressure offset in one-dimensional compressible model

To illustrate evolution of the mean pressure inside the computational domain, we consider 1D compressible dry model with a localized heat source Q mimicking latent heating in the moist multidimensional case. The governing set of equations has the form

 
formula
 
formula
 
formula
 
formula

where ρ(x, t), u(x, t), T(x, t), and p(x, t) are 1D fields of the air density, velocity, temperature, and pressure, respectively; Rd is the gas constant of the dry air; d/dt ≡ ∂/∂t + u∂/∂x; and ≡ ∂/∂x. Model variables are split into a constant reference state and perturbations:

 
formula
 
formula
 
formula
 
formula

and the set of Eqs. (B1)(B4) is solved for the perturbations with all of them initially set to zero. The reference state is defined as ρ0 = 1.16 kg m−3, u0 = 0 m s−1, T0 = 300 K, and p0 = 105 Pa. The domain length is 20 km with the grid length of Δx = 31.25 m. The prescribed heat source is Q = 0.1 K s−1 in the 2-km-wide region in the middle of the domain and zero otherwise. Model equations are converted into the flux form and solved applying the nonoscillatory forward-in-time method using the 1D MPDATA advection scheme (cf. Smolarkiewicz and Szmelter 2009). The density and temperature equations are solved first and the equation of state is used to obtain pressure needed to solve the momentum equation. Model equations are integrated for 240 s applying 0.12-s time step.

Two types of lateral boundary conditions are considered: “open” and periodic. The open boundary scheme applies an absorber defined in Davies (1976) over six points near the lateral boundaries to remove perturbations before they reach the boundaries; reference-state values are prescribed at the boundaries. Such an approach mimics open boundaries where perturbations are allowed to propagate freely out of the domain. In contrast, periodic lateral boundary conditions result in a closed domain because perturbations are never able to escape: as they leave one side of the domain, they come back from the other.

Numerical solutions to the 1D heating problem show strong increase (decrease) of the temperature (density) in the central region of the domain where the prescribed heating is applied, a weak flow away from the heating (maximum around 0.3 m s−1), and solitary sound waves propagating away from the central region and allowing adjustments of model fields away from the prescribed heating. The adjustment through propagating solitary waves appears similar to a case of the heating in motionless stratified atmosphere (e.g., Nicholls et al. 1991), although obviously solitary sound waves rather than gravity bores occur in our case. The solitary sound waves are either absorbed near lateral boundaries when open boundary conditions are used, or they return through the other boundary in the case of periodic conditions. The two situations result in a different mean pressure in the domain.

Figure B1 compares solutions at the end of the simulations for the two types of boundary conditions. The figure shows spatial distributions of relative perturbations of the temperature (T/T0), density (ρ/ρ0), pressure (p/p0), and the horizontal velocity—the latter scaled by the speed of sound cs = (cp/cυ × p0/ρ0)1/2 = 347 m s−1 and multiplied by 20 to fit the vertical scale of the figure. In both cases, the temperature increase within the heating region leads to the reduction of the density, with the pressure field remaining approximately constant across the domain. The key difference between the two solutions is the magnitude of the pressure field. In the case of periodic conditions, the pressure increases from the initial value of 1000 hPa to about 1011 hPa, consistent with the increase of the mean temperature (i.e., 〈p′〉 = 〈ρRdT′〉, where the angle brackets depict domain average). When open boundary scheme is used, the pressure field increases just slightly (about 0.8 hPa), illustrating that the scheme only approximately mimics unobstructed propagation of perturbations out of the domain. The domain-averaged density does not change when periodic conditions are used (i.e., the total mass is conserved), but it does change when open boundary scheme is used. The latter is because of the mass exchange through the lateral boundaries.

Fig. B1.

One-dimensional compressible model solutions to the problem of local heating for open and periodic lateral boundary conditions.

Fig. B1.

One-dimensional compressible model solutions to the problem of local heating for open and periodic lateral boundary conditions.

Overall, the simple 1D test highlights the role of lateral boundary conditions in the formation of the mean pressure in the domain for the compressible system of equations. This aspect is relevant to the comparison among anelastic and compressible solutions as discussed in the main text.

b. Hydrostatic adjustment in anelastic and compressible models

We consider anelastic and compressible solutions to a problem of horizontally uniform heating in the two-dimensional nonrotating stratified dry atmosphere. Since the atmosphere and the heating are horizontally uniform, the problem is one dimensional and, in principle, the model described in the previous section can be used after adding stratification in Eq. (B5) to Eq. (B8) and gravitational acceleration in Eq. (B2). However, we use the full EULAG model in a 2D configuration to retain numerical implementation of both anelastic and compressible versions. The atmosphere is assumed to have a constant static stability of 1.3 × 10−5 m−1. The domain size is 2 × 10 km2 (width × height) with a uniform grid length of 100 m. The heat source is Q = 5 K h−1 in a layer between 1 and 3 km and zero otherwise. Surface temperature and pressure are 1000 hPa and 300 K, respectively. Periodic lateral boundary conditions and rigid/free-slip boundary conditions at the top and bottom of the domain are applied. The models are integrated for 1 h with 0.25-s time step.

The atmosphere remains motionless in the anelastic simulation. The compressible model simulates insignificant vertical motions of the atmosphere (velocities ~10−2 m s−1) resulting from the hydrostatic adjustment involving vertically propagating sound waves. The anelastic and compressible models start from a balanced (i.e., satisfying compatibility conditions) initial state with no pressure perturbations. Without any heat sources (i.e., Q = 0) the initial state remains unchanged, and no pressure perturbations are generated. Hence, the only change in pressure comes from the nonzero temperature source mimicking the effects of the latent heating in moist simulations.

Evolution of the pressure perturbation profiles for the anelastic and compressible simulations are shown in Fig. B2. Because the nonhydrostatic pressure perturbations are not allowed to develop in the simulations (except for those minuscule due to sound waves in the compressible system), the only change of the pressure field comes from a gradual buildup of the hydrostatic component in response to the temperature increase. Slight changes of the temperature and density profiles away from the prescribed heating are simulated in the compressible model (not shown). These come from adjustments of the temperature and density profiles outside the heating region facilitated by vertically propagating sound waves. For periodic boundary conditions, conservation of mass in the compressible model leads to a small increase of the density outside the heating region to compensate the density decrease within the heating region. The density profile in the anelastic model does not change by definition.B1 To allow a better comparison of the anelastic and compressible solutions, the pressure offset (Fig. B2c) that develops similarly as in 1D tests was subtracted from the compressible solution.

Fig. B2.

Profiles of horizontally averaged pressure perturbations for (a) ANES and (b) COMP models for the problem of the horizontally uniform heating applied between 1 and 3 km. (c) Time evolution of the pressure offset for the compressible model (black) that was subtracted from the COMP profiles in (b). In addition, theoretical estimation of the offset based on the assumption that the pressure increases owing to gradual temperature increase within a closed-box domain is shown.

Fig. B2.

Profiles of horizontally averaged pressure perturbations for (a) ANES and (b) COMP models for the problem of the horizontally uniform heating applied between 1 and 3 km. (c) Time evolution of the pressure offset for the compressible model (black) that was subtracted from the COMP profiles in (b). In addition, theoretical estimation of the offset based on the assumption that the pressure increases owing to gradual temperature increase within a closed-box domain is shown.

The pressure offset simulated in the compressible model can also be estimated theoretically for both anelastic and compressible models. We use the same formulation of the equation of state as used in the model; that is,

 
formula

where p0 = 105 Pa. For the compressible model, the offset is derived by averaging local values of the pressure derived from Eq. (B9)—that is, including spatial variability of the density and potential temperature fields. For the anelastic model, only spatial variability of the potential temperature is considered and the reference-state density is used. Nevertheless, the estimates are accurate as shown in Fig. B2c. It should be stressed that the development of the pressure offset cannot be simulated in the anelastic system because the elliptic solver adds no constant to the pressure perturbation solution (see appendix in KSG13). The fact that the offset can be recovered for the anelastic model has potentially important implications.

To verify the origin of pressure perturbations shown in Fig. B2, a hydrostatic component of pressure perturbations p′ can be calculated for each model based on its hydrostatic balance equation. The compressible hydrostatic balance is

 
formula

where ρ is the air density, ρ′ = ρρ0 denotes the density perturbation around the initial profile ρ0, and g is the gravitational acceleration. The anelastic hydrostatic balance is

 
formula

where and define reference (and initial) profiles of the density and potential temperature, respectively, and θ′ is the potential temperature perturbation about .

Integrating Eq. (B10) upward or downward and applying p′ at the starting level derived from Eq. (B9) as p′ = pp0 leads to pressure profiles exactly as shown in Fig. B2b for the compressible model. Similarly, integrating Eq. (B11) upward and downward from the level where p′ is zero in the anelastic system (around 2.4 km) results in the pressure perturbation profiles as shown in Fig. B2a. Interestingly, if the information about an actual pressure perturbation from the compressible model is used to calculate the anelastic hydrostatic balance [see Eq. (B11)], then the hydrostatic profile from the compressible model (i.e., Fig. B2b) is approximately recovered. Therefore, the key differences between the anelastic and compressible hydrostatic adjustments originate from different boundary conditions for the pressure field rather than slightly different formulations of the hydrostatic balance.

In summary, the differences in the pressure field simulated by anelastic and compressible models can be understood through the development of different hydrostatic components in both model formulations. These components have negligible impact on moist processes because of its small amplitude, at least for the cases considered in the current study.

REFERENCES

REFERENCES
Bannon
,
P. R.
,
1996
:
On the anelastic approximation for a compressible atmosphere
.
J. Atmos. Sci.
,
53
,
3618
3628
, doi:.
Bryan
,
G. H.
, and
J. M.
Fritsch
,
2002
:
A benchmark simulation for moist nonhydrostatic numerical models
.
Mon. Wea. Rev.
,
130
,
2917
2928
, doi:.
Davies
,
H. C.
,
1976
:
A lateral boundary formulation for multilevel prediction models
.
Quart. J. Roy. Meteor. Soc.
,
102
,
405
418
, doi:.
Durran
,
D. R.
,
1989
:
Improving the anelastic approximation
.
J. Atmos. Sci.
,
46
,
1453
1461
, doi:.
Grabowski
,
W. W.
,
1998
:
Toward cloud resolving modeling of large-scale tropical circulations: A simple cloud microphysics parameterization
.
J. Atmos. Sci.
,
55
,
3283
3298
, doi:.
Grabowski
,
W. W.
, and
P. K.
Smolarkiewicz
,
1990
:
Monotone finite-difference approximations to the advection-condensation problem
.
Mon. Wea. Rev.
,
118
,
2082
2098
, doi:.
Grabowski
,
W. W.
, and
P. K.
Smolarkiewicz
,
1996
:
On two-time-level semi-Lagrangian modeling of precipitating clouds
.
Mon. Wea. Rev.
,
124
,
487
497
, doi:.
Grabowski
,
W. W.
, and
P. K.
Smolarkiewicz
,
2002
:
A multiscale anelastic model for meteorological research
.
Mon. Wea. Rev.
,
130
,
939
956
, doi:.
Held
,
I. M.
, and
M. J.
Suarez
,
1994
:
A proposal for the intercomparison of the dynamical cores of the atmospheric general circulation models
.
Bull. Amer. Meteor. Soc.
,
75
,
1825
1830
, doi:.
Klemp
,
J. B.
, and
R. B.
Wilhelmson
,
1978
:
The simulation of three-dimensional convective storm dynamics
.
J. Atmos. Sci.
,
35
,
1070
1096
, doi:.
Kurowski
,
M. J.
,
B.
Rosa
, and
M. Z.
Ziemianski
,
2011
:
Testing the anelastic nonhydrostatic model EULAG as a prospective dynamical core of a numerical weather prediction model. Part II: Simulations of supercell
.
Acta Geophys.
,
59
,
1267
1293
, doi:.
Kurowski
,
M. J.
,
W. W.
Grabowski
, and
P. K.
Smolarkiewicz
,
2013
:
Toward multiscale simulation of moist flows with soundproof equations
.
J. Atmos. Sci.
,
70
,
3995
4011
, doi:.
Lipps
,
F. B.
, and
R. S.
Hemler
,
1982
:
A scale analysis of deep moist convection and some related numerical calculations
.
J. Atmos. Sci.
,
39
,
2192
2210
, doi:.
Margolin
,
L. G.
,
P. K.
Smolarkiewicz
, and
Z.
Sorbjan
,
1999
:
Large-eddy simulations of convective boundary layers using nonoscillatory differencing
.
Physica D
,
133
,
390
397
, doi:.
Miglietta
,
M. M.
, and
R.
Rotunno
,
2005
:
Simulations of moist nearly neutral flow over a ridge
.
J. Atmos. Sci.
,
62
,
1410
1427
, doi:.
Nicholls
,
M. E.
,
R. A.
Pielke
, and
W. R.
Cotton
,
1991
:
Thermally forced gravity waves in an atmosphere at rest
.
J. Atmos. Sci.
,
48
,
1869
1884
, doi:.
Prusa
,
J. M.
,
P. K.
Smolarkiewicz
, and
A. A.
Wyszogrodzki
,
2008
:
EULAG, a computational model for multiscale flows
.
Comput. Fluids
,
37
,
1193
1207
, doi:.
Smolarkiewicz
,
P. K.
,
2006
:
Multidimensional positive definite advection transport algorithm: An overview
.
Int. J. Numer. Methods Fluids
,
50
,
1123
1144
, doi:.
Smolarkiewicz
,
P. K.
, and
L. G.
Margolin
,
1998
:
MPDATA: A finite-difference solver for geophysical flows
.
J. Comput. Phys.
,
140
,
459
480
, doi:.
Smolarkiewicz
,
P. K.
, and
A.
Dörnbrack
,
2008
:
Conservative integrals of adiabatic Durran’s equations
.
Int. J. Numer. Methods Fluids
,
56
,
1513
1519
, doi:.
Smolarkiewicz
,
P. K.
, and
J.
Szmelter
,
2009
:
Iterated upwind schemes for gas dynamics
.
J. Comput. Phys.
,
228
,
33
54
, doi:.
Smolarkiewicz
,
P. K.
,
C.
Kühnlein
, and
N.
Wedi
,
2014
:
A consistent framework for discrete integrations of soundproof and compressible PDEs of atmospheric dynamics
.
J. Comput. Phys.
,
263
,
185
205
, doi:.
Temam
,
R.
,
2006
:
Suitable initial conditions
.
J. Comput. Phys.
,
218
,
443
450
, doi:.
Weisman
,
M. L.
, and
J. B.
Klemp
,
1982
:
The dependence of numerically simulated convective storms on vertical wind shear and buoyancy
.
Mon. Wea. Rev.
,
110
,
504
520
, doi:.
Wilhelmson
,
R.
, and
Y.
Ogura
,
1972
:
The pressure perturbation and the numerical modeling of a cloud
.
J. Atmos. Sci.
,
29
,
1295
1307
, doi:.

Footnotes

*

The National Center for Atmospheric Research is sponsored by the National Science Foundation.

1

Strictly speaking, such a statement is correct for the anelastic model only. Using the anelastic continuity equation, one can show that the vertical gradient of the domain-averaged vertical velocity has to vanish. Since the vertical velocity is zero at the surface, the average has to vanish at all levels. For the compressible system with periodic lateral boundary conditions, horizontally averaged vertical velocity can only be nonzero because of the flow compressibility.

2

This result has been substantiated by a complimentary pseudoincompressible experiment that closely reproduces the compressible solution (cf. section 4.1 in SKW14). Unlike the anelastic system, a pseudoincompressible system (Durran 1989) is based on the unapproximated evolutionary form of the momentum equation.

A1

To facilitate the implementation of the framework for dry and moist equations, the moist compressible solver assumes θd/θdaθ/θa while formulating the buoyancy force; typically neglecting O(10−4) departures from unity—since the density potential temperature θd = θ(1 + qυ/ϵ)/(1 + qt), and θda is its ambient value.

B1

As shown in Bannon (1996), one can reconstruct the thermodynamic density within the heating region for the anelastic system using a linearized diagnostic relation between density, pressure, and the potential temperature perturbations [see Eq. (3.17) therein]. However, the reconstruction cannot account for the density increase outside the heating region.