## Abstract

Moist anelastic and compressible numerical solutions to the planetary baroclinic instability and climate benchmarks are compared. The solutions are obtained by applying a consistent numerical framework for discrete integrations of the various nonhydrostatic flow equations. Moist extension of the baroclinic instability benchmark is formulated as an analog of the dry case. Flow patterns, surface vertical vorticity and pressure, total kinetic energy, power spectra, and total amount of condensed water are analyzed. The climate benchmark extends the baroclinic instability study by addressing long-term statistics of an idealized planetary equilibrium and associated meridional transports. Short-term deterministic anelastic and compressible solutions differ significantly. In particular, anelastic baroclinic eddies propagate faster and develop slower owing to, respectively, modified dispersion relation and abbreviated baroclinic vorticity production. These eddies also carry less kinetic energy, and the onset of their rapid growth occurs later than for the compressible solutions. The observed differences between the two solutions are sensitive to initial conditions as they diminish for large-amplitude excitations of the instability. In particular, on the climatic time scales, the anelastic and compressible solutions evince similar zonally averaged flow patterns with the matching meridional transports of entropy, momentum, and moisture.

## 1. Introduction

This is the third paper in a series of works devoted to investigating the relative merits of the anelastic and compressible moist dynamics across the range of scales, from small to meso to planetary. The first paper (Kurowski et al. 2013) introduced our notion of an all-scale moist simulation and documented the consistent anelastic and compressible solutions for idealized shallow convective and orographic cloud formations. The second paper (Kurowski et al. 2014) analyzed moist deep convection and demonstrated that anelastic approximation accurately represents severe convective dynamics. The current paper extends the earlier studies to planetary scales. Normal-mode analysis (Davies et al. 2003) indicates that the anelastic approximation of Lipps and Hemler (1982) misrepresents large-scale atmospheric flows compared to predictions based on fully compressible Euler equations. On the other hand, multiple-scales asymptotic analyses (Dolaptchiev and Klein 2009, 2013) show that, at synoptic–planetary length and time scales, atmospheric motions are predominantly anelastic. Both results are correct, and they do not contradict each other. Depending on the focus of interests, perturbations about a predominantly anelastic state of the atmosphere can be judged sufficiently large to disprove the suitability of an anelastic model for, for example, weather prediction and climate studies (Davies et al. 2003). However, practical suitability limits—or alternatively, manifestations of the unsuitability in simulations of realistic atmospheric flows—have not been established, especially for moist global flows.

For weather and climate simulations, sound waves are energetically unimportant but can be computationally demanding. Consequently, the soundproof systems of PDEs that retain thermal aspects of compressibility but analytically filter out rapidly propagating sound waves have certain appeal for nonhydrostatic modeling. The advancement of high-performance computing over the last two decades has enabled development of large-scale, high-resolution models. This, in turn, revived discussions about the range of validity of soundproof approximations (e.g., Cullen et al. 2000; Davies et al. 2003; Klein et al. 2010), originally proposed as an alternative to the compressible Euler equations for limited-area applications (cf. Lipps and Hemler 1982; Durran 1989). Notwithstanding the low impact of sound waves on atmospheric circulations, there are important issues associated with the use of soundproof equations in large-scale modeling. These are outlined below to set the ground for the subsequent discussion of the anelastic and compressible solutions at the planetary scale.

The soundproof systems are built on linearizations discarding certain perturbational terms that are arguably small. Historically, scale analysis arguments employed in derivations of soundproof approximations limited their validity to weak background stratifications, thus questioning the utility of soundproof systems for simulation of realistic atmospheres even on meso-*γ* scales. Recently, Klein et al. (2010) showed the formal validity of the anelastic (Lipps and Hemler 1982) and pseudoincompressible (Durran 1989) systems for realistic background stratifications with potential temperature variations of 30–50 K across the troposphere. Accordingly, Kurowski et al. (2014) demonstrated close agreement between anelastic and compressible numerical solutions for moist deep convection, with differences between results for two different mathematical formulations insignificant compared to sensitivities to numerical details and subgrid-scale parameterizations. Similarly, the results of Smolarkiewicz et al. (2014) for orographically induced stratospheric gravity waves also showed excellent agreement between compressible and soundproof solutions. As far as the largest horizontal scales are concerned, the assumption of a horizontally homogeneous base state, inherent in the anelastic system, yields maximum meridional temperature deviations on the order of 20% when compared to the midlatitude profiles (e.g., Held and Suarez 1994, hereafter HS94). This is perceived as significant, although such deviations are still an order of magnitude smaller than the reference values. More importantly, linearization of the pressure gradient term in the evolutionary form of the anelastic momentum equation abbreviates baroclinic production of vorticity. This has been demonstrated in Smolarkiewicz et al. (2014), where related departures of the anelastic solutions from the pseudoincompressible and compressible results were shown to be significant for synoptic scales.

A normal-mode analysis (Davies et al. 2003; Arakawa and Konor 2009) for linearized equations on the Cartesian *f* plane reveals key differences between compressible and soundproof dispersion relations, especially for the longest internal atmospheric modes. Short and mesoscale horizontal modes (up to ~100 km) are correctly represented by the soundproof approximations. The differences begin to appear for long horizontal (~1000 km) and deep vertical (~40 km) modes. The linear theory predicts energy redistribution between the modes for the pseudoincompressible approach and a positive phase shift for the anelastic approach (Davies et al. 2003). Smolarkiewicz et al. (2014) illustrated some of these predictions in a baroclinic instability simulation for an idealized global atmosphere. Notwithstanding a good agreement between pseudoincompressible and compressible solutions, they also showed that soundproof systems yield higher group velocity: the pseudoincompressible wave propagated about 0.5 m s^{−1} faster than the compressible one, with a similar difference between the anelastic and pseudoincompressible solutions.

Soundproof systems dictate the solution of elliptic Poisson equations for pressure perturbations about a balanced hydrostatic ambient state to ensure mass continuity of the resulting flow. These perturbations are typically assumed to be small and excluded from the model thermodynamics (cf. appendix A in Lipps and Hemler 1982). In principle, moist processes such as saturation adjustment are affected by this simplification. A heuristic analysis of Kurowski et al. (2013) shows that the small-scale nonhydrostatic component of the pressure perturbations is typically less important than the larger-scale quasi-hydrostatic component. Furthermore, their results demonstrate that the anelastic nonhydrostatic pressure perturbations compare well with the compressible counterparts and thus are suitable for reconstructing the full pressure field. Nonetheless, numerical experiments with small-scale cloud dynamics and orographic flows (Kurowski et al. 2013) and with moist deep convection (Kurowski et al. 2014) documented the negligible impact of the pressure perturbations on moist thermodynamics.

In the current study, two planetary-scale dry benchmarks are extended to take into account effects of moist processes, including phase changes and precipitation. The two benchmarks simulate, respectively, the formation of a baroclinic wave following Jablonowski and Williamson (2006, hereafter JW06) and the idealized climate of HS94. The extended JW06 problem epitomizes the development of midlatitude weather systems that involve intense vorticity dynamics and carry large amounts of moisture and energy over distances of thousands of kilometers. The calculations compare short-term (a week or so) deterministic solutions that develop from a localized smooth perturbation. In contrast, the moist HS94 problem addresses long-term statistical properties of the equilibrium climate and accompanying meridional transports.

## 2. The consistent soundproof/compressible numerical framework

The research tool employed in the study is the all-scale Eulerian/semi-Lagrangian (EULAG) model designed to integrate four different dynamical-core equation sets—the anelastic (Lipps and Hemler 1982), the pseudoincompressible (Durran 1989), and two distinct adaptations of compressible Euler equations (to be specified shortly)—with minimal differences in the numerics (Smolarkiewicz et al. 2014). Here, we continue in the spirit of the earlier works (Kurowski et al. 2013, 2014) and focus on anelastic and compressible simulations. For inviscid adiabatic problems, the pseudoincompressible system is a relatively straightforward extension of the anelastic equations (Smolarkiewicz and Dörnbrack 2008) and gives results close to compressible Euler equations for a broad range of scales (Smolarkiewicz et al. 2014). Generally, however, numerical integrations of the pseudoincompressible equations are more involved because of the distinctive elliptic constraint that explicitly includes the diabatic source (Almgren 2000; O’Neill and Klein 2014; Duarte et al. 2015). The global problems addressed in the current paper feature multiplicity of intermittent localized heat sources associated with evolving cloud fields. This complicates the integrability of the pseudoincompressible elliptic constraint and deserves a separate study. Nevertheless, the pseudoincompressible integrations with diabatic contributions neglected in the elliptic constraint [cf. section 6 in Durran (1989) for a discussion] closely match the compressible results, in the spirit of adiabatic simulations in Smolarkiewicz et al. (2014).

The numerical design for the soundproof and compressible mathematical formulations of the governing PDEs in EULAG (Smolarkiewicz et al. 2014; Kurowski et al. 2014) provides a consistent framework operating on the same set of dependent variables, written in essentially the same perturbation form, and using the same two-time-level Eulerian/Lagrangian principal integration algorithm (viz. time stepping) for all soundproof and compressible dynamical cores. In all formulations, the principal algorithm shares the same advection scheme (for all prognostic variables) and the same elliptic solver. Furthermore, all dynamical cores share the same curvilinear coordinate transformations, computational grid, spatial discretization, and parallelization schemes.

All prognostic equations are cast into the conservative flux form and integrated using the nonoscillatory forward-in-time approach. In the default model algorithm, the rotational, buoyant and acoustic modes are all treated implicitly,^{1} admitting the same large time steps in compressible and soundproof systems. In the acoustic variant of the compressible model, the thermodynamic pressure is diagnosed directly from the potential temperature and density, as opposed to the elliptic boundary value problem employed in the large time step models; whereas the rotational and gravitational modes are still treated implicitly as in the large time step models (Smolarkiewicz et al. 2014). The key difference between the soundproof and the implicit compressible models is in the form of the elliptic boundary value problem. The implicit compressible model solves the Helmholtz problem composed of three Poisson operators akin to the operator employed in the adiabatic soundproof systems (Smolarkiewicz et al. 2014). The extension of the implicit compressible model to incorporate effects of moisture (including contributions to the Helmholtz elliptic equation) has recently been discussed in Kurowski et al. (2014). The three systems employed in this study are the compressible implicit (COMP), compressible explicit (COMPe; i.e., acoustic) and the anelastic (ANES).

Because the study focuses on the comparison of dynamical cores rather than on a comprehensive representation of moist processes, only large-scale condensation/evaporation is considered, and neither convection nor boundary layer parameterizations are used. To keep the setup simple and easy to reproduce, the Kessler warm rain parameterization is employed with the autoconversion threshold of 0.5 g kg^{−1}. The bulk moist thermodynamics is the same as in Kurowski et al. (2014). Ice forming processes are excluded from the model setup. The compressible model employs full pressure in the saturation adjustment unless otherwise stated, whereas the soundproof moist thermodynamics is based on either the ambient pressure or the full pressure, with the latter including pressure perturbations as in the generalized anelastic model of Kurowski et al. (2013, 2014). For the use in moist thermodynamics and/or analysis, the soundproof pressure perturbation needs to be filtered out from the unphysical component related to the null space of the discrete nabla operator employed in the momentum equation.

## 3. Baroclinic wave development

### a. Simulation setup

The original JW06 setup assumes a steady-state atmosphere of the Earth-like rotating sphere with two midlatitude zonal jets symmetrical across the equatorial *x*–*z* plane. In the jets, the maximum zonal velocity is 35 m s^{−1}, and the prescribed flow is in geostrophic and thermal balance. Baroclinic instability is triggered by an isolated solenoidal perturbation of the zonal flow for the Northern Hemisphere jet, with the center of the perturbation located at 40°N, 20°E.

The horizontal mesh consists of 256 × 128 grid points on the regular longitude–latitude grid. The vertical domain of 23 km is covered with 48 uniformly distributed levels and the corresponding grid interval m. Rayleigh damping is applied in the vicinity of the poles to suppress development of superresolved modes resulting from the convergence of meridians. No explicit diffusion is used. The time step for COMP and ANES is 300 s, and it is 2 s for COMPe. Selected sensitivity experiments with COMP also employ the acoustic time step. The base state for the soundproof models has a constant static stability of 1.02 × 10^{−5} m^{−1}, close to the tropospheric value in the midlatitudes.

In this paper, we consider two alternative ways for extending an established dry benchmark to the moist atmosphere, each with merits on its own. The first alternative that adds moisture as a deviation to a known dry setup (Grabowski and Smolarkiewicz 2002; Park et al. 2013) will be used for the HS94 idealized climate in section 4. In this section, the second alternative is used that accounts for the moist phase (Waite and Snyder 2013; Kurowski et al. 2014) while adjusting the dry setup so as to maintain the same ambient balance in the moist case. Here, the adjustment maintains the geostrophic and thermal balance of the dry benchmark important for controlled growth of baroclinic instability. In effect, the moist extension of the dry JW06 setup retains the original initial fields of zonal wind and pressure, while the potential temperature field *θ* and the density *ρ* are altered. In particular, *θ* is adjusted such that the density potential temperature,

is equal to the initial potential temperature of the original dry case (, , and are the water vapor, cloud water, and rainwater mixing ratios, respectively; , with *R* and denoting the gas constants for the dry air and water vapor). Assuming initial , a known , and a constant surface pressure equal to the reference value hPa, and integrating upward the hydrostatic balance relation (metric terms aside),

assures that the moist Exner pressure is equal to the initial Exner pressure of the original dry case ( is specific heat of dry air at constant pressure, and *g* gravitational acceleration). In turn, the geostrophic balance

assures the equality of the zonal flow velocities of the dry and moist flow, , and thus the thermal wind balance

of the dry setup ( and correspond to differentiation in the north–south and vertical directions, and *f* is the Coriolis parameter). Furthermore, to satisfy the gas law in the moist case, the initial dry density is evaluated directly from

where , and *θ* is already adjusted as explained above. The adjusted potential temperature is about 2% colder near the surface at the equator, with the percentage decreasing both poleward and with height. The adjusted dry density behaves similarly, with *ρ* also about 2% lower near the surface at the equator. Notably, the results are independent of the initial amount of moisture, as long as phase changes are turned off.

To prescribe a realistic initial field of , we start with a zonally and meridionally homogeneous relative humidity RH field, prescribed constant in the lower troposphere and smoothly transitioning to zero aloft

where is the value of RH below the height of , and is the level above which there is no water vapor. The values = 2 km, = 6 km, and were selected for simulations discussed in this paper.^{2} Having defined relative humidity in (6) allows us to prescribe . For simplicity, to avoid an iterative adjustment of *θ*, is calculated from RH assuming the original *θ* from the dry setup, which, in turn, results in a slightly modified RH field compared to (6) (Waite and Snyder 2013).

To examine the effects of moisture on the large-scale flows and to assess the role of pressure perturbations in moist thermodynamics, numerical experiments with baroclinic waves are conducted with increasing levels of complexity. Starting with the same initial conditions, three different setups are considered: (i) phase changes switched off (setup C1); (ii) only condensation/evaporation allowed (setup C2); and (iii) both condensation and precipitation processes included (setup C3).

### b. Flow patterns

When phase changes of the water substance are excluded, liquid water remains zero throughout the simulation, and water vapor becomes a passive scalar. In such a situation, potential temperature has formally no sinks or sources. This configuration enables verification of the numerical model, since dry and moist solutions should develop in the same way, assuming that the effects of moisture, other than the 2% adjustments of *θ* and *ρ* discussed above, are negligible in this specific case. Figure 1 shows the surface virtual potential temperature, , and pressure perturbations at day 8 for COMP, COMPe, and ANES simulations. The differences between moist anelastic and compressible models reflect those for (explicitly dry) simulations in Fig. 2 of Smolarkiewicz et al. (2014). In particular, the baroclinic wave propagates significantly faster, while the entire wave train develops significantly more slowly, in the ANES model.

Figure 2 complements Fig. 1 with zonal cross sections at 53°N through the surface fields of the virtual potential temperature and pressure perturbations and the water vapor mixing ratio. The faster zonal wave propagation in the ANES model is apparent. A detailed inspection shows that the phase shift is the smallest on the left-hand side of the wave train (i.e., in the tail of the wave structure). The initial perturbation develops into the most mature eddy, and each subsequent wave that develops in the ANES model has a smaller phase shift. At day 8, the differences of the wave packet leading edge location (i.e., at about 220°E) correspond to about 1–1.5 m s^{−1} faster propagation of the ANES solution as compared to the compressible result. Consistent with Fig. 1, the magnitude of pressure perturbations is uniformly lower as the structure develops more slowly. Comparing magnitudes of and perturbations at this zonal cross section is inconclusive, as it is obscured by the frontogenesis seen in the troughs of the pressure wave in Fig. 1 for both compressible solutions.

Including condensation/evaporation (setup C2) has a small impact on patterns of surface virtual potential temperature and pressure (not shown). Similarly, adding precipitation (setup C3) has negligible impact on the surface flow patterns, as illustrated in Fig. 3. For the small- and mesoscale dynamics, precipitation is a driving factor in the formation of cold pools and squall lines, and it affects the large-scale flow. In our case, however, the flow is driven by synoptic-scale wave dynamics, and both condensation and precipitation have only a minor impact on the flow. Other fields, such as the surface pressure perturbations and vorticity, are practically the same for the three setups C1, C2, and C3. The latter is substantiated with Fig. 4 displaying vertical (relative) vorticity and flow vectors superimposed with isolines of pressure perturbations for the setup C3. The patterns of the virtual potential temperature perturbations (not shown) are essentially the same as those for the potential temperature perturbations in the dry case. In other words, given the initial balance [(2)–(4)], the solution for the surface buoyancy perturbations does not depend on humidity of the air. In fact, the dry and moist EULAG solutions for a range of are practically undistinguishable (not shown). The presence of condensation and precipitation does not seem to affect the phase shift documented in Fig. 2.

A vertical cross section at 53°N through the fields of vertical velocity and cloud water mixing ratio for the most complete C3 setup (i.e., with phase changes and precipitation included), is shown in Fig. 5. Vertical motions are about a few centimeters per second, and they span the entire depth of the troposphere. Waves tilt westward—that is, in the direction opposite to the flow (line A in the top panel of Fig. 5)—which is typical for horizontally propagating baroclinic waves (Holton 1979, chapter 6). Clouds form on fronts separating cold and warm air masses that tilt eastward (line B in the top panel of Fig. 5). The frontal cloud structure in this region is about 7 km deep for COMP and COMPe. The anelastic model yields similar patterns but with vertical currents and cloud fields significantly underdeveloped. Note that all these structures evince many small-scale details that are typically smoothed by horizontal diffusion or semi-implicit schemes with large time steps (cf. Figs. 8 and 12 in Polvani et al. 2004). In particular—as documented by animations of the COMP solutions using time steps 2, 50, 150, and 300 s (not shown)—the small-scale features apparent in the upper troposphere for the COMPe solution resolve the gravity wave response to localized heat release evolving at the time scale similar to the Brunt–Väisälä period ( s).

### c. Surface vorticity and the growth of eddies

One key difference between the anelastic and compressible PDEs lies in the momentum equation. The (perturbation) pressure gradient force is linearized and potential in the anelastic approximation of Lipps (1990). Consequently, the horizontal gradient of buoyancy is the sole baroclinic source of vorticity, directly producing (breeze like) circulations in vertical planes. In compressible Euler equations, the form of the pressure gradient force implies the nonlinear baroclinic source of vorticity that can directly produce circulations in horizontal planes. In small-scale dynamics, buoyant vorticity production dominates, rendering the baroclinic production of vertical vorticity negligible. However, this is not necessarily the case for planetary scales, as shown in Smolarkiewicz et al. (2014) for the dry atmosphere. Moisture adds another dimension to this argument (Cao and Cho 1995), because the baroclinic source also accounts for water vapor and liquid water mixing ratios included in and [cf. (A1)–(A6) in Kurowski et al. 2014].

To evaluate the role of the nonlinear vorticity production in evolution of the baroclinic instability, the history of the maximum surface vertical vorticity for COMP, COMPe, and ANES is displayed in Fig. 6a. All three dynamical cores evince roughly the same slow vorticity growth in the first five days of the evolution. Notably, in this first stage of the evolution (referred to as “linear,” following Prusa and Gutowski 2010) the anelastic results match closely the original JW06 integrations of the hydrostatic primitive equations (Smolarkiewicz 2011). After day 5, the vorticity growth in compressible solutions suddenly accelerates, marking the onset of the nonlinear, rapid-growth phase for the fastest growing baroclinic eddy. This transition is controlled by the baroclinic source of vorticity, as substantiated with the ad hoc experiments of Smolarkiewicz et al. (2014; see the last paragraph of their section 4.1) demonstrating reproducibility of compressible (and pseudoincompressible) results by arbitrarily manipulating the coefficients of the pressure gradient force. Between the days 8 and 9, the fastest growing eddy reaches its maximum strength, the wave breaks, and the regular structure of the wave train begins to disintegrate. The flow transitions into a strongly nonlinear multiscale regime characteristic of geophysical turbulence. Further analysis of the maximal becomes ambiguous, as it does not identify anymore the distinct eddy. Figure 6b complements the history of maximal (Fig. 6a) with the history of maximal surface meridional velocity, sometimes used to illustrate the baroclinic wave evolution. In contrast to the vorticity, it does not discriminate between the two phases of the evolution before the wave breaking, showing a steady exponential growth since day 1 (cf. Fig. 4 in Park et al. 2013).

An alternative view on the instability is presented in Fig. 7. The figure shows the history of the selected surface virtual potential temperature isolines around the day of transition to the nonlinear phase (day 5.2) for COMP (setup C3). Meridional perturbations of are relatively small and regular before day 5.2. When baroclinic vortices start to affect the large-scale flow, the stirring process gradually intensifies, and so does meridional displacement of the surface temperature isolines.

As for the maximum , the main difference between ANES and COMP/COMPe seems to be in the starting time of the rapid-growth phase. The anelastic model needs roughly one additional day to form the large-scale perturbations from which baroclinic eddies further develop. Once the rapid growth sets in at day 6, the history of the maximum closely resembles that for COMP and COMPe. However, the anelastic growth rate is 15%–20% smaller, and the rapid-growth phase lasts for about 3.8 days: that is, roughly half a day longer than for COMP/COMPe. The difference develops mostly during the last two days of the rapid-growth phase, with the growth rate in the first two days almost the same for all three models. The wave breaking in the anelastic model occurs for about a 15% lower value of . Arguably, all these differences are because of different baroclinic vorticity production in ANES (Smolarkiewicz et al. 2014). A small enhancement of for COMPe around the wave breaking is the effect of using a smaller time step, as the COMP solution with the acoustic time step evinces virtually the same behavior (not shown).

Latent heat release increases the growth rate and shortens by about 3–4 h the duration of the rapid-growth phase. This mechanism of moist invigoration of synoptic systems has already been discussed by several authors (e.g., Gutowski et al. 1992; Reed et al. 1993; Booth et al. 2013). Although the history of the maximum documents this effect, detecting it from the flow patterns in Fig. 3 is hardly possible. The solutions with and without precipitation are almost the same (not shown).

Auxiliary simulations conducted on a coarser grid of 2.8° × 2.8° revealed that the duration of the linear and rapid-growth phases of the wave evolution and the rate of the rapid growth are sensitive to integration accuracy. In particular, coarser integrations resulted in the spin-up time of 4 (COMP) or 4.8 (ANES) days, as opposed to 5.2 or 6 days, whereas the rapid-growth phase lasted for about 4.2 (COMP) or 5 days (ANES), as opposed to 3.4 or 3.8 days, respectively. Simulations of complex atmospheric flows applying various governing equations are typically performed with various numerical models and solution methods, thus leaving the origin of discrepancies uncertain (cf. section 7b in Ullrich et al. 2014, for a discussion). The use of a single numerical model with consistent numerics in all simulations bolsters our confidence in the integrity of the results. The most significant difference noted so far, between the anelastic and compressible solutions, is in the duration of the linear phase. It is not unreasonable, therefore, to anticipate that in realistic meteorological situations, where cyclogenesis starts from relatively large perturbations, the disparity between the solutions will manifest differently. We shall return to this point later in the paper.

### d. Eddy kinetic energy and minimum surface pressure

A standard way of evaluating the strength of baroclinic eddies is by analyzing the evolution of eddy kinetic energy (EKE) and minimum surface pressure (e.g., Lorenz 1955; Simmons and Hoskins 1978; Pavan et al. 1999; Booth et al. 2013; Ullrich et al. 2014). Here, EKE measures the magnitude of the flow perturbations with respect to the ambient state, whereas the minimum pressure together with the maximum mark the center of the most developed eddy, at least in the evolution’s linear phase.

The history of the global EKE integral is shown in Fig. 8a. Although diluted by global integration, the three stages of the evolution corresponding to the history of in Fig. 6 still can be identified. In this metric, however, the first stage lasts about 10 h longer in ANES than in COMP and COMPe. In contrast to the maximum (Fig. 6a), the transition from the slow spinup to the rapid growth is gradual, and it takes about a day or so. The most evident disparity between ANES and COMP/COMPe develops during the rapid-growth phase, which has two distinct growth rates for the anelastic and compressible systems. This stage ends when the wave breaking sets in. The EKE increases more rapidly until day 8 for COMP and COMPe and until day 10 for ANES. At later times, the growth rates are almost the same for all models, and the differences between the anelastic and compressible solutions developed at the previous stage remain at the same level. Note that the EKE increases over the entire 14-day integration. This is the main difference between global simulations and a typical baroclinic life cycle simulation in a periodic channel (cf. Figs. 3a,b in Pavan et al. 1999). The latter offers a more controlled environment in which a baroclinic eddy reaches its maximum strength within a few days and then gradually decays. In our experiments, however, a positive tendency lasts for about a month (not shown) as perturbations spread across the entire planet, with the onset of circulations in the Southern Hemisphere in the third week (not shown). Only then the EKE begins to decrease. The influence of moisture on EKE is relatively small, as are the regions of cloud presence. For the rapid-growth phase, latent heat release enhances EKE by a few percent at most.

The histories of the minimum surface pressures are depicted in Fig. 8b. The ANES pressure perturbation starts to decrease about one day later than those of COMP and COMPe, with the latter two almost perfectly matching each other. The presence of condensation invigorates the dynamics of baroclinic eddies as the surface pressure minimum drops by an additional 5–10 hPa.

### e. Kinetic energy spectra

Observations (e.g., Nastrom and Gage 1985; Lindborg 1999) and numerical experiments (e.g., Hamilton et al. 2008; Skamarock et al. 2014) suggest that the canonical energy spectrum for well-developed atmospheric circulations is proportional to at synoptic scales and shallows to at the mesoscale (here, *k* denotes the horizontal wavenumber). The comparison of tropospheric kinetic energy spectra for COMP, COMPe, and ANES simulations at days 10 and 14 is depicted in Fig. 9. The spectra are derived from the horizontal wind perturbations and averaged over the lowest 8 km of the troposphere, excluding surface values. In our simulations, it takes several days for the atmosphere to develop the steady-state energy cascade across the entire range of scales. The longest modes accumulate most of the energy, which is then transferred through the synoptic scales down to the mesoscale, starting from approximately wavenumber 10. COMP, COMPe, and ANES spectra remain consistent at the largest scales (i.e., for wavenumbers from 1 to 5). Significantly less energy is accumulated in the ANES synoptic-to-mesoscale cascade owing to a slower development of baroclinic eddies.

At day 10, the three spectra already display some features of the canonical spectrum. Because the ANES solution lags behind, the steady-state energy cascade has not yet been developed, and at mesoscale the spectrum tends to follow the −3 slope, rather than the canonical . At day 14, the solutions have accumulated more energy in smaller scales, and their characteristics are closer to the canonical spectrum. Moreover, the shape of the ANES spectra closely resembles that for COMP. For all three models, the slope of the synoptic part is slightly steeper than −3, consistent with the results of Skamarock et al. (2014). Similarly to the total eddy kinetic energy, the amount of energy for the ANES cascade at day 14 is close to that for COMP and COMPe at day 10. Note also that the COMPe model features more energy at the finest scales than COMP, arguably because of a better temporal resolution of small-scale features with the acoustic time step.

Simulated moisture effects enhance the kinetic energy mainly at smaller scales (). Additionally, the buoyancy production due to latent heat release helps to establish multiscale flow sooner, as it earlier reshapes the tails of power spectra toward the slope (cf. Fig. 9a). Generally, these tendencies agree with the Waite and Snyder (2013) results for midlatitude *f*-plane simulations in a rectangular periodic channel. Here, the moist spectra are also shallower than the slope at the smallest scales, where, in the absence of an explicit subgrid-scale closure, the spectral behavior is controlled by the model nonoscillatory numerics [Domaradzki et al. 2003; see also Schaefer-Rolffs and Becker (2013) for a related discussion]. The overall history of the spectra indicates that the synoptic–mesoscale break is associated with the transition from two- to three-dimensional turbulence as the spectra start evolving from a slightly greater than −3 slope for both synoptic- and mesoscale ranges, and the mesoscale range gradually evolves toward the slope after the onset of the wave-breaking phase.

### f. Condensed water and the role of pressure perturbations

Figure 10 shows the history of the total cloud water amount for the COMP, COMPe, and ANES simulations with and without precipitation included in the calculations. Condensation first occurs during the fifth (COMP and COMPe) or sixth (ANES) day of simulations, although it seems to start later for the C2 setup as Figs. 10a and 10b use differently scaled ordinates. For the case with precipitation (C3), only a small fraction of the condensed water is accumulated in clouds, and most of it is converted into rain and falls out quickly. For the case without precipitation (C2), the only mechanism of reducing liquid water is evaporation, which is less efficient than rain formation and fallout.

When the rain autoconversion is excluded, all solutions evolve smoothly, and the water amount keeps increasing throughout the simulation. The COMP and COMPe yield almost the same cloud water amount, with only minor differences in the last day of the simulations (i.e., for the turbulent phase). The ANES solution grows significantly more slowly, with the difference increasing with time. At day 14, the ANES cloud amount reaches about one-third of that for COMP and COMPe. Accounting for pressure perturbations in moist thermodynamics has only a minor impact on the solutions.

Adding precipitation impacts the integrations in several ways. First, the simulations become more variant, especially in the turbulent phase. Second, the simulations are more sensitive to the choice of a time step, as the maximal differences between COMP and COMPe can reach about 10% and, on average, COMPe evinces somewhat lower cloud water amount (i.e., precipitates more efficiently). Third, including pressure perturbations in moist thermodynamics impairs the relatively regular history simulated in the C2 setup, making the results less conclusive. The sensitivity to the choice of a time step turns out to be significantly larger than to including/excluding pressure perturbations. As compared to the C2 case, the differences between the compressible and anelastic model formulations somewhat diminish, and the ANES model features about 50% of the total water amount compared to the COMP solution at day 14. Similarly to the total eddy kinetic energy, the ANES values at day 14 are close to those for COMP and COMPe at day 10.

### g. Realistic initial state

The experimental setup of JW06 assumes that the initial zonally symmetric flow is balanced and laminar, and the instability grows from a small Gaussian perturbation. Such an idealization of the initial state makes the JW06 setup easy to implement in a broad range of numerical models. However, it obscures theoretical/numerical model comparisons by elevating the role of a particular initial condition that enables growth of the solution disparities already in the linear spin-up phase. Although this emphasis on the slow incubation of the most unstable modes is important in itself, we also wish to assess differences of anelastic and compressible dynamics in applications akin to developed weather. Consequently, we conduct additional simulations based on a more realistic initial state. To eliminate differences arising in the spin-up phase, we use the compressible solution with the C3 setup at day 5.2 (marking the onset of the rapid-growth phase; Fig. 6) as the initial condition for ANES.

Figure 11 shows the history of the surface maximum vertical vorticity for ANES initiated as described above, together with the COMP solution already presented in Fig. 6. The growth rate of the maximum for ANES is larger than in Fig. 6, and the wave breaking takes place a few hours after COMP. Although the maximum for ANES is 15% lower than for COMP at the onset of the wave breaking, it does not imply that the metric remains uniformly smaller, as documented for later integration times.

Figure 12 presents the surface virtual potential temperature and pressure perturbations for ANES, as well as their zonal cross sections along 53°N at day 8. Although less pronounced, the overall disparity of the anelastic solution and the corresponding COMP results in Figs. 1 and 3 is still apparent. Zonal cross sections document the average difference in the wave amplitude reaching 20%, and the phase shift between COMP and ANES reduced in proportion to the integration time. The total EKE and minimum surface pressure in COMP and ANES show closer agreement (Fig. 13), yet the difference between ANES and COMP grows in time but at a smaller rate than in Fig. 8. For the surface pressure the trend reverses at later times, and the pressure is uniformly lower in ANES than in COMP after day 12, already in the turbulent phase.

## 4. Moist idealized climate

### a. Simulation setup

The HS94 climate benchmark was proposed to facilitate intercomparison of dynamical cores of general circulation models. The simulated global circulation is driven by the Newtonian relaxation of the temperature field to the prescribed zonally symmetric radiative equilibrium and the near-surface Rayleigh damping mimicking the boundary layer friction.

Following Smolarkiewicz et al. (2001), the numerical setup assumes a 32-km-deep domain resolved with 41 levels. Vertical stretching mimics a uniform grid in pressure coordinates corresponding to a fixed exponential profile with a 7-km height scale and a 25-hPa increment. In the horizontal, the surface of the sphere is resolved with 128 × 64 uniform latitude × longitude grid points. All simulations span 1000 days, resolved with a 120-s uniform temporal increment for dry and moist experiments. The base-state atmosphere assumes constant static stability of 1.02 × 10^{−5} m^{−1}. The absorbing sponge layer extends from *z* = 15.4 km to the model top. The inverse time scale of the absorber increases linearly from zero at the bottom of the layer to 1 day at the top. A small explicit diffusion of 30 m^{2} s^{−1} is applied in momentum equations. Results from the first 200 days are considered a spin-up time and excluded from analysis. Moist extension of the benchmark is in the spirit of Grabowski and Smolarkiewicz (2002), with moisture added as a deviation to the original dry setup. Here, the forcing term for the water vapor mixing ratio follows Newtonian relaxation of the temperature field, with the equilibrium determined from (6) using the equilibrium *θ* of the dry setup. For simplicity, only warm rain processes are considered. COMP and ANES models use either the base-state pressure profile in moist thermodynamics or include pressure perturbations to reconstruct full pressure, as described below.

### b. Preamble: Dry solutions

Figure 14 displays zonally averaged climatological means of the potential temperature, zonal wind, EKE, and meridional transports of the entropy and zonal momentum for the dry experiment. Here, all perturbation variables are evaluated with respect to climatological means. Both models simulate consistent potential temperature distribution and zonal jet structures, with slightly stronger equatorial stratification in the upper troposphere for ANES. The maximum strength of jets in the ANES models is a few percent smaller than in COMP (see Table 1 for maxima of various quantities in the dry and moist HS94 simulations). Small differences can also be detected for the westerlies and subpolar circulations, which are both weaker in ANES. The patterns of meridional transports and EKE are similar for both models. Furthermore, all the transports are also in quantitative agreement between the two solutions.

### c. Climatological means and transports

Zonally averaged means from moist simulations are depicted in Fig. 15. Because the tropics are now well mixed for moist saturated air, a stronger stratification of the potential temperature is simulated. As in the dry case, fields of the potential temperature and zonal jets are similar between the two models. Although a weaker polar circulation can still be detected for ANES, many of the differences documented for the dry setup decrease in the moist simulation. Unlike in the baroclinic instability study of section 3, the equilibrium potential temperature is the same for the dry and moist setups. Consequently, the water vapor amplifies meridional gradients of the density potential temperature. This, in turn, roughly doubles the kinetic energy of the synoptic–planetary-scale modes (not shown), upon which the moist models simulate about 1.4 times stronger midlatitude jets (cf. Table 1). Furthermore, meridional transports get notably stronger, though their general patterns compare well with the dry counterparts. The enhancement of meridional transports is most likely caused by stronger baroclinicity of the atmosphere and thus more energetic midlatitude synoptic-scale eddies, which are primary agents of the poleward advection of heat, momentum, and moisture (e.g., Arakawa 1975). Moreover, the kinetic energy for moist solutions features an additional maximum in the tropics around the top of the troposphere. The ANES results remain in good agreement with the COMP solutions, with similar maxima for most of the zonal means as shown in Table 1.

### d. Distribution and transport of moisture

Zonally averaged distribution and transport of atmospheric water vapor and mean surface precipitation are presented in Fig. 16. The climatological equilibrium moisture distribution is virtually the same for ANES and COMP. The fluxes of the water vapor mixing ratio, as well as the meridional transports of entropy and zonal momentum, are also in good agreement, with similar structures around the equatorial belt and in the subtropical-to-midlatitude zones. The figure also presents climatological means for pressure perturbations, which indicate the quasi-hydrostatic component that builds up in a response to persistent low-tropospheric heating around the tropics. Once the global circulation is fully developed after about 60 days of spinup, this distribution only minimally changes in time. The quasi-hydrostatic component is the main contribution to the pressure perturbation, as the local changes due to the development of baroclinic eddies are typically an order of magnitude smaller and have both positive and negative excursions. Nevertheless, including the perturbations in moist thermodynamics results in only minimal changes in the surface precipitation (cf. Fig. 16). More detailed effects of the pressure perturbation on the climatological means are documented in Table 1. Although zonal jets remain practically unchanged, the global maxima for meridional transports and kinetic energy are generally several percent larger when the full pressure is used in moist thermodynamics. On the other hand, cumulative surface precipitation—a sink term for the atmospheric moisture—does not seem to be sensitive to the choice of the governing equations and whether pressure perturbations are included in moist thermodynamics. One needs to keep in mind that the moisture forcing is nonconservative, because its magnitude is proportional to the difference between the actual state and the prescribed equilibrium. In Earth’s atmosphere, however, the moisture budget is controlled by the partitioning of the surface water exchange between evaporation and precipitation (Arakawa 1975; Peixoto and Oort 1984).

Figure 17 shows Hovmöller diagrams of the surface precipitation from the tropics (i.e., at the equator) and from the midlatitudes (at 55°N). COMP and ANES solutions show similar surface precipitation patterns. In midlatitudes, the precipitation zones are associated with frontal systems and are driven by an eastward propagation of baroclinic eddies. In the tropics, the surface precipitation pattern features more complex organization, with large-scale (wavenumbers 4 and 5) coherent structures propagating westward and embedded strongly precipitating deep convective systems. The degree to which the equatorial convection is organized varies with time, and the large-scale pattern can become more chaotic after a period of well-organized propagation. The strength of the midlatitude precipitation seems to oscillate with time as well.

Normalized PDFs of the surface precipitation for the tropics and midlatitudes are depicted in Fig. 18. They are first calculated for eight 100-day time windows (such as that from Fig. 17) and then averaged over an 800-day period of time, with a standard deviation based on the 100-day means. The PDFs closely agree between the simulations. For the midlatitudes (frontal precipitation) the differences between the models remain smaller than one standard deviation, whereas for the tropics (deep convection), they are hardly perceptible.

## 5. Concluding remarks

Scale- and normal-mode analyses of the Euler equations for the dry atmosphere indicate that soundproof approximations are well suited for representing small- to mesoscale atmospheric flows. Beyond that, it is generally agreed that predictive capabilities of soundproof approximations diminish. However, the actual magnitude of the solution disparities accounting for uncertainty of simulated natural phenomena has never been conclusively demonstrated for nonlinear flows. This is especially true for the moist precipitating atmosphere, the theoretical analyses of which are rare and usually limited to small scales. This paper investigates the relative performance of the anelastic approximation for synoptic- and planetary-scale moist flows. The consistent numerical framework for integrating the soundproof and compressible equations of Smolarkiewicz et al. (2014) and Kurowski et al. (2014) is used to compare the anelastic solutions against the corresponding compressible results. The dry benchmarks of planetary baroclinic instability and idealized climate are extended to account for moist processes, with an aim to provide minimal models for natural weather and climate. The baroclinic wave benchmark aids to quantify the relative capability of numerically integrated anelastic and compressible PDEs for deterministic forecasts. The climate benchmark extends the baroclinic wave study to the equilibrium climate and addresses the relative capability of the two PDE systems for predicting mean flows statistics and the associated meridional transports.

Numerical solutions to the baroclinic wave evolution demonstrate an important role of the baroclinic vorticity production. In the anelastic system, the pressure-perturbation gradient force is linearized and potential, and the anelastic model simulates slower development of baroclinic instability. This contrasts with both the pseudoincompressible (not shown) and fully compressible solutions, the governing equations for which include the unapproximated pressure gradient force. Furthermore, the numerical results indicate that the initial growth rate of surface vorticity—a measure of the strength of synoptic-scale eddies—is similar in all systems considered. The main difference is a delayed onset of the rapid-growth phase in the anelastic model. During the rapid-growth phase, the anelastic growth rate is 15% smaller than its compressible counterpart, and the wave breaking occurs at 15% lower values of the maximum surface vorticity. It is noteworthy that, when the anelastic model uses a realistic initial state with large-amplitude perturbations, the differences between the two models are less apparent.

For the anelastic and compressible equations, the coarsely resolved climate equilibria compare well between dry and moist setups. In the moist case, zonal jets and meridional transports reach higher maxima, but a general picture and the conclusions remain intact. Meridional distribution of entropy and moisture, as well as the structure of zonal jets, are all in good agreement between compressible and anelastic models. Only small differences appear for meridional transports of the entropy, zonal momentum, and moisture. The impact of the dynamic pressure and density perturbations on moist solutions was found to be small as well. The climate results highlight a potential issue related to the nonconservative design of the Held–Suarez benchmark, since meridional means seem unaffected by the differences in poleward transports. Because of that, more tests are needed to clarify the suitability of soundproof approximation for global climate simulations, perhaps applying aquaplanets. Furthermore, desirable are similar comparisons at much higher resolutions (hardly affordable to the authors at the present time) that may enhance the role of baroclinic vorticity production in moist simulations. If this is the case, then statistics targeting extreme events are also necessary to better quantify the role of moisture and the degree of anelasticity at climatic scales.

An advantage of the EULAG framework is its numerical consistency. In the course of simulations, we have found a large sensitivity of the integrations to the numerical aspects, including accuracy of the advection scheme, time-step size, resolution, numerical filters, and tuning parameters. It should be stressed, however, that the differences documented in the paper remain valid for all configurations tested as long as the numerical solvers for integrating different equations consequently employ the same set of controlling parameters in the numerical framework. In contrast, the different numerical environments can easily obscure solutions’ disparities because of inherent differences between the theoretical model formulations. This is highlighted in Fig. 19, which juxtaposes the C1 setup solutions from Fig. 6 with the corresponding compressible result that employs heavy filtering in the spirit of the composite schemes of Liska and Wendroff (1998).^{3} In the maximum surface vorticity metric, the filtered compressible solution closely matches the anelastic result.

## Acknowledgments

Personal reviews from Drs. Christian Kühnlein and Nils Wedi are gratefully acknowledged. Comments from Prof. Rupert Klein and an anonymous referee helped to improve the presentation. This material is based upon work partially supported by the U.S. Department of Energy under Award DE-SC0006748, and it was partially done during MJK’s stay at NCAR. Part of the research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. 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).

## REFERENCES

*Workshop on developments of the ECMWF model*, Reading, United Kingdom, ECMWF. [Available online at http://nldr.library.ucar.edu/repository/collections/OSGC-000-000-010-248.]

*An Introduction to Dynamic Meteorology.*Academic Press, 391 pp.

*Proc. Fifth European Conf. on Computational Fluid Dynamics (ECCOMAS CFD 2010)*, Lisbon, Portugal, European Community on Computational Methods in Applied Sciences, 1453. [Available online at http://congress.cimne.com/eccomas/cfd2010/papers/01453.pdf.]

*Proc. ECMWF Workshop on Nonhydrostatic Modelling*, Reading, UK, ECMWF, 15 pp. [Available online at http://old.ecmwf.int/publications/library/ecpublications/_pdf/workshop/2010/Nonhydrostatic_modelling/Smolarkiewicz.pdf.]

## Footnotes

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

^{1}

For large-scale problems, fast propagating gravity waves severely limit the model time step if solved explicitly (Smolarkiewicz et al. 2001; Grabowski and Smolarkiewicz 2002; Smolarkiewicz 2011). For instance, explicit anelastic integrations of the baroclinic instability problem addressed in this paper require 20-times-smaller time steps than the corresponding implicit integrations; see Fig. 4 and the accompanying discussion in Smolarkiewicz (2011). Furthermore, the consistent trapezoidal-rule time integration of all principal forcings also enhances model accuracy (Dörnbrack et al. 2005; Wedi and Smolarkiewicz 2006).

^{2}

Larger values of lead to large-scale convective overturning in the tropics, because the potential instability in the lower troposphere (viz. the equivalent potential temperature decreasing with height) makes the moist layer convectively unstable when brought to saturation.

^{3}

Namely, the nonoscillatory advection algorithm of the EULAG model (Smolarkiewicz et al. 2014) was set to use the generic first-order-accurate upwind scheme every third time step.