## Abstract

Many high-resolution atmospheric models can reproduce the qualitative shape of the atmospheric kinetic energy spectrum, which has a power-law slope of −3 at large horizontal scales that shallows to approximately −5/3 in the mesoscale. This paper investigates the possible dependence of model energy spectra on the vertical grid resolution. Idealized simulations forced by relaxation to a baroclinically unstable jet are performed for a wide range of vertical grid spacings Δ*z*. Energy spectra are converged for Δ*z * 200 m but are very sensitive to resolution with 500 m ≤ Δ*z* ≤ 2 km. The nature of this sensitivity depends on the vertical mixing scheme. With no vertical mixing or with weak, stability-dependent mixing, the mesoscale spectra are artificially amplified by low resolution: they are shallower and extend to larger scales than in the converged simulations. By contrast, vertical hyperviscosity with fixed grid-scale damping rate has the opposite effect: underresolved spectra are spuriously steepened. High-resolution spectra are converged except for the stability-dependent mixing case, which are damped by excessive mixing due to enhanced shear over a wide range of horizontal scales. It is shown that converged spectra require resolution of all vertical scales associated with the resolved horizontal structures: these include quasigeostrophic scales for large-scale motions with small Rossby number and the buoyancy scale for small-scale motions at large Rossby number. It is speculated that some model energy spectra may be contaminated by low vertical resolution, and it is recommended that vertical-resolution sensitivity tests always be performed.

## 1. Introduction

The horizontal wavenumber spectrum of kinetic energy in the atmosphere has a distinct double power-law shape. The spectral slope of −3 at large scales shallows to approximately −5/3 below scales of a few hundred kilometers (Gage 1979; Nastrom and Gage 1985; Cho et al. 1999). Over the last two decades, numerical models have become increasingly good at reproducing the general shape of this spectrum; these include global models (e.g., Koshyk and Hamilton 2001; Hamilton et al. 2008; Brune and Becker 2013; Skamarock et al. 2014) and regional models (e.g., Skamarock 2004; Waite and Snyder 2009, 2013; Bei and Zhang 2014; Peng and Zhang 2014; Peng et al. 2015). The inability of some models to capture the shallowing of the mesoscale spectrum has led to insights about numerical damping and subgrid parameterization (Shutts 2005; Augier and Lindborg 2013). As a result, the kinetic energy spectrum has become a useful diagnostic for assessing the ability of models to capture the cascade of kinetic energy through the mesoscale.

The resolution of the shallow mesoscale part of the spectrum requires fine horizontal meshes. In computational studies of the energy spectrum, emphasis has therefore been placed on decreasing the horizontal grid spacing Δ*x*. Global models with Δ*x* ~ *O*(1) km are now possible, which allow for a full decade of the mesoscale spectrum to be resolved (Skamarock et al. 2014). On the other hand, the issue of vertical resolution has not received as much attention. A wide range of vertical grid spacings Δ*z* are used in the studies cited above. For example, in the upper troposphere, Δ*z* ≈ 1500 m in Hamilton et al. (2008) and ≈100 m in Waite and Snyder (2013). Furthermore, most studies of model energy spectra have not analyzed the sensitivity of results to vertical resolution. An exception is Brune and Becker (2013), who compared spectra obtained with typical low vertical resolution (30 levels, with Δ*z* ≈ 1500 m in the upper troposphere) and much higher resolution (100 levels, with Δ*z* ≈ 200 m). Surprisingly, they found the shallow mesoscale spectrum in the upper troposphere to be more pronounced in the low-vertical-resolution case.

The question of what is appropriate vertical resolution in atmospheric simulations has been addressed in many studies. For quasigeostrophic (QG) motions away from boundaries, the grid should be consistent with the QG aspect ratio *f*/*N*, where *f* is the Coriolis parameter and *N* is the buoyancy frequency (Charney 1949). As a result, the vertical grid spacing must be fine enough to resolve horizontal-grid-scale QG vortices (e.g., Lindzen and Fox-Rabinovitz 1989); that is,

Finer vertical resolution may be necessary in the presence of fronts (Pecnick and Keyser 1989; Persson and Warner 1991; Snyder et al. 1993; Iga et al. 2007), gravity wave critical layers (Lindzen and Fox-Rabinovitz 1989), or to resolve surface temperature advection by QG vortices (Tulloch and Smith 2009; Bembenek et al. 2015). These studies suggest that inadequate vertical resolution may lead to spurious gravity waves or noise at small horizontal scales.

The condition in (1) assumes QG dynamics. However, for sufficiently small Δ*x* well inside the mesoscale, grid-scale motion will not be QG, and (1) may not be appropriate. Sub-QG motions may include stratified turbulence, gravity waves, and moist convection, which have different requirements for vertical resolution. Stratified turbulence, which has a −5/3 energy spectrum like that observed in the mesoscale (e.g., Lindborg 2006; Riley and Lindborg 2008), develops shear layers with thickness around the buoyancy scale *L*_{b} ≡ 2*πU*/*N*, where *U* is a representative horizontal velocity scale (Billant and Chomaz 2001; Waite and Bartello 2004). Typical mesoscale values of *L*_{b} are from a few hundred meters to 1 km when the 2*π* factor is included (e.g., Brune and Becker 2013). The numerical grid must resolve these layers, leading to a different condition on vertical resolution:

Simulations of rotating–stratified turbulence at different *f* and *N* show that the characteristic vertical scale makes a transition from the QG scale to the buoyancy scale as the Rossby number Ro exceeds *O*(1) (Lindborg 2005; Waite and Bartello 2006). This transition suggests that, as Δ*x* is reduced, the condition on Δ*z* may change from (1) to (2) once the grid-scale Ro exceeds an *O*(1) threshold.

In addition to the number of levels, the treatment of damping and mixing at small vertical scales varies between models, and the influence of different approaches on the horizontal spectrum has not been explored. The Weather Research and Forecasting (WRF) Model (Skamarock et al. 2008) is typically run with an upwind biased advection scheme, which introduces some numerical damping at small vertical scales (Skamarock 2004). For example, a third-order scheme introduces a fourth-order numerical hyperviscosity (∂^{4}/∂*z*^{4}) with a coefficient proportional to the local Courant number (Wicker and Skamarock 2002). Other models use energy-conserving vertical differencing schemes [e.g., the AGCM for the Earth Simulator (AFES) model; Ohfuchi et al. 2004]. In addition, explicit vertical mixing may be applied. In the free troposphere, the vertical eddy viscosity and diffusivity are usually computed locally and depend on the Richardson number Ri. This dependence suppresses or eliminates vertical mixing when conditions are stable (e.g., Mellor and Yamada 1982; Hong et al. 2006).

In idealized simulations of homogeneous stratified and rotating–stratified turbulence, hyperviscosity with fixed horizontal and vertical coefficients are more common than flow-dependent eddy viscosity (e.g., Waite and Bartello 2006; Lindborg 2006; Kitamura and Matsuda 2006; Bartello 2010; Deusebio et al. 2013). In such cases, the coefficients are chosen to be as small as possible while providing sufficient grid-scale dissipation to avoid numerical instabilities and spectral blocking. Idealized simulations often have much higher vertical resolution than the global and regional model simulations cited above; for example, Deusebio et al. (2013) used 1024 levels in the vertical. Stratified turbulence simulations without rotation that respect the vertical-resolution condition in (2) have a horizontal wavenumber spectral slope close to −5/3 (Lindborg 2006). However, when (2) is violated, vertical hyperviscosity damps large horizontal scales, inhibiting the energy cascade and resulting in a very steep energy spectrum (e.g., Waite and Bartello 2004).

In this paper, we investigate the dependence of horizontal wavenumber energy spectra on vertical resolution. Idealized simulations of rotating–stratified turbulence, forced by a baroclinically unstable jet, are presented. The influence of different vertical mixing schemes is considered, including no explicit mixing, hyperviscosity, and stability-dependent eddy viscosity. This work is motivated by an apparent contradiction in previous work: stratified turbulence, hypothesized to be a model for the mesoscale −5/3 spectrum, requires resolution of the buoyancy scale in the vertical; on the other hand, atmospheric models are routinely able to capture a shallow mesoscale spectrum with much lower vertical resolution. The work presented here shows that underresolution in the vertical can profoundly affect the horizontal spectrum, and in some cases even amplify it, but the details depend on the treatment of vertical mixing. Our approach is outlined in section 2, results are presented in section 3, and conclusions are given in section 4.

## 2. Method

We consider an idealized dry dynamical core governed by the Boussinesq equations on an *f* plane with constant buoyancy frequency:

Here **u** = (*u*, *υ*, *w*) is the velocity; ∇ and ∇_{h} are the three-dimensional and horizontal gradient operators, respectively; *p* is the pressure scaled by a reference density; *b* is the buoyancy, which is proportional to the fluctuation potential temperature; **F**_{u} is the velocity forcing; **D**_{u} and *D*_{b} are the vertical mixing terms for velocity and buoyancy, respectively; and is the unit upward normal. Periodic boundary conditions are employed in the horizontal and vertical directions (note, however, that only the potential temperature perturbation, and not the full potential temperature, is periodic in *z*). The domain size is *L* × *L* × *H*, where *L* = 4000 km and *H* = 20 km. Unless otherwise stated, we take *f =* 10^{−4} s^{−1} and *N* = 10^{−2} s^{−1}.

Simulations are forced by relaxation to zonal jets with a sinusoidal meridional and vertical structure:

where

The maximum velocity is *U* = 30 m s^{−1}, and the dimensionless jet wavenumbers are *k* = 2 and *m* = 1. Here and throughout, we refer to dimensionless (integer) wavenumbers, which are made dimensional by multiplication by the horizontal and vertical wavenumber spacings Δ*k* = 2*π*/*L* and Δ*m* = 2*π*/*H*. The jet size, from maximum eastward to maximum westward velocity, is *L*/4 = 1000 km in the horizontal plane and *H*/2 = 10 km in the vertical direction. The Rossby and Froude numbers of the reference jet are therefore

The relaxation time scale is *τ* = 1 day.

The forcing in (6) is applied to all wavevectors **k** = (*k*_{x}, *k*_{y}, *k*_{z}) inside the cube |*k*_{x}|, |*k*_{y}|, and |*k*_{z}| ≤ 2, so large-scale modes that are not part of the jet are relaxed to zero. This large-scale damping prevents accumulation of energy at the largest scales. The forced jet is baroclinically unstable, so synoptic and mesoscale features develop naturally via growth and saturation of this instability. While idealized, this yields an arguably more realistic forcing than random excitation of Fourier modes (e.g., Waite and Bartello 2006; Deusebio et al. 2013).

Horizontal and vertical mixing are treated separately. In the horizontal plane, fourth-order hyperviscosity is employed with constant coefficients *ν*_{h} = *κ*_{h} for momentum and buoyancy, respectively. The value of *ν*_{h} is chosen to be as small as possible while preventing the build-up of grid-scale energy. The grid-scale damping rate is around 22 min. In the vertical direction, three possibilities for the momentum and buoyancy mixing **D**_{u} and *D*_{b} are considered. In the first case, no vertical mixing is applied. In the second case, hyperviscosity with fixed coefficients is employed; that is,

where again *ν*_{υ} = *κ*_{υ}. The coefficients are varied with resolution to keep the grid-scale damping rate the same as in the horizontal. Last, vertical Smagorinsky (1963) mixing modified by a simple stability function is considered:

where the eddy viscosity and diffusivity are given by

ℓ = 30 m is the mixing length,

is the vertical shear,

is the local Richardson number, and

is the stability function (following Holtslag and Boville 1993; Neale et al. 2010). Note that *F*(Ri) grows unbounded for Ri → −∞ but approaches zero rapidly for Ri > 0 (Fig. 1). Therefore, unstable motions are strongly damped, but stable motions are only weakly affected.

Equations (3)–(5) are solved with a Fourier-based spectral code. Since spectral discretization introduces no numerical viscosity, all diffusion is controlled by the horizontal hyperviscosity and vertical mixing schemes. Fast Fourier transforms are used to evaluate nonlinear terms on a physical-space grid, and aliasing errors are eliminated with the two-thirds rule (Orszag 1971). The number of physical-space grid points is *n*_{h} × *n*_{h} × *n*_{υ}. After dealiasing, the dimensionless horizontal and vertical truncation wavenumbers are

rounded to the nearest integer. The effective grid spacings are therefore

Time stepping is third-order Adams–Bashforth, and hyperviscosity terms are treated implicitly with a Crank–Nicolson approach. The vertical eddy viscosity, when present, is integrated explicitly.

Horizontal resolution is fixed in all cases at *n*_{h} = 512, which corresponds to *k*_{T} = 170 and Δ*x* = 11.8 km. The horizontal hyperviscosity coefficient is correspondingly fixed at *ν*_{h} = 2.5 ×10^{25} m^{8} s^{−1}. A range of vertical resolutions from *n*_{υ} = 16 (Δ*z* = 2.0 km) to *n*_{υ} = 256 (Δ*z* = 118 m) are considered. Simulations with *n*_{υ} = 16 and *n*_{υ} = 32 are spun up from initial conditions at *t* = 0 given by the reference jet in (7)–(8) plus low-level random noise. Simulations with higher vertical resolution are restarted from the *n*_{υ} = 32 case at *t* = 12 days. Simulations with vertical eddy viscosity in (11) are restarted from simulations with no vertical mixing at *t* = 12 days. All simulations are run to *t* = 24 days, and energy spectra and other diagnostics are averaged over the last 9 days of integration.

Table 1 summarizes the parameters and scalar quantities for all cases. An overview of the simulations with no vertical mixing is given in Fig. 2, which shows time series of average kinetic energy per unit mass. In the simulations started from *t* = 0, the baroclinic instability grows over several days and saturates around *t* = 8 days. High-resolution simulations with *n*_{υ} ≥ 64 diverge from one another a few days after being restarted at *t* = 12 days.

Kinetic energy spectra are computed from the three-dimensional discrete Fourier transforms of the velocity field. The kinetic energy per unit mass associated with dimensionless wavevector **k** is

where the caret denotes the Fourier coefficient and the asterisk denotes the complex conjugate. The horizontal and vertical integer wavenumbers {*k*_{x}, *k*_{y}} and *k*_{z} are between ±*k*_{T} and ±*m*_{T}, respectively. One-dimensional horizontal wavenumber spectra are constructed by summing over all *k*_{z} and over annuli of unit width on the *k*_{x}–*k*_{y} plane (e.g., Waite and Bartello 2004); that is,

where *k*_{h} ≥ 0 are integers and Summing over *k*_{z} is equivalent to computing *k*_{h} spectra at every level and averaging in the vertical, which is often done in more realistic models. The binning over integer wavenumbers in (20) leads to bumps in the spectra due to sampling issues at the annulus interfaces. This is corrected by multiplying (20) by the ratio of the annulus area to the actual number of modes in each annulus, which smoothes the *k*_{h} spectra slightly.

## 3. Results

### a. No vertical mixing

Figure 3a shows kinetic energy spectra for the runs with no vertical mixing at different vertical resolutions. All spectra have roughly the same qualitative shape as the Nastrom and Gage (1985) spectrum: they are peaked at wavenumbers 1 and 2 associated with the large-scale forcing, the peak is followed by a power law with a slope around −3, and the slope shallows significantly at larger mesoscale wavenumbers. Beyond *k*_{h} ≈ 100, at length scales of a few grid points, the spectrum falls off rapidly due to the horizontal hyperviscosity. The large-scale spectral amplitude varies slightly with vertical resolution, due to the differences in average kinetic energy over the averaging interval (Fig. 2). This variance is corrected for in Fig. 3b by normalizing the spectra with the energy in *k*_{h} = 3, which allows for a clearer comparison of the mesoscale spectra in each case.

Despite the broadly similar structure of the spectra in Fig. 3, the shallowing of the spectrum at mesoscale wavenumbers is actually highly sensitive to vertical resolution. A wavenumber characterizing the transition from a −3 to a shallower slope is identified by finding the local minimum of (after smoothing the spectra with a Gaussian filter with width *σ* = 2). These transition wavenumbers, which we call *k*_{m}, are given in Table 1. With *n*_{υ} = 16, the spectrum shallows at *k*_{m} = 14, which corresponds to a wavelength of 286 km. Downscale of this transition, a broad mesoscale power law is visible for 20 80. The slope of the mesoscale spectrum in this case is approximately −0.8, which is notably shallower than −5/3. However, as *n*_{υ} increases, the transition wavenumber increases and the corresponding amplitude of the mesoscale spectrum decreases. As a result, the shallow mesoscale spectrum is broader and more pronounced at low vertical resolution [consistent with the results in Brune and Becker (2013)]. The spectra are reasonably converged at our highest resolutions, although there are some minor differences as *n*_{υ} increases from 128 to 256. Interestingly, even in the converged simulations, the spectral slope of the shallow part of the spectrum is much shallower than −5/3. A power-law fit over 60 ≤ *k*_{h} ≤ 100 in the *n*_{υ} = 256 case gives a slope of −0.8.

To understand the dependence of the transition wavenumber *k*_{m} on vertical resolution, consider the QG resolution condition in (1). According to this criterion, for simulations that are under-resolved in the vertical, the maximum dimensional horizontal wavenumber for which QG motions can be consistently captured is *f*/*N* times the maximum vertical wavenumber. Expressed in terms of dimensionless wavenumbers, this gives a maximum horizontal wavenumber of

the values of which are given in Table 1. For the two lowest vertical resolutions *n*_{υ} = 16 and 32, the transition wavenumber *k*_{m} is close to the maximum QG wavenumber in (21): 14 versus 10 for *n*_{υ} = 16, and 18 versus 20 for *n*_{υ} = 32. For higher vertical resolutions, while the QG wavenumber in (21) continues to increase for increasing *n*_{υ}, the transition wavenumber *k*_{m} appears to converge to a value around 25, corresponding to a wavelength of 160 km.

Convergence of the spectrum requires that the vertical grid be fine enough to resolve the vertical scales associated with the full range of horizontal scales. To identify these scales, first consider the spectral Rossby number associated with a particular horizontal wavenumber *k*_{h},

which is shown in Fig. 4a for the high-resolution case with *n*_{υ} = 256. At small and intermediate *k*_{h} with Ro(*k*_{h}) < 1, the dynamics are expected to be approximately QG, in which case the vertical wavenumber associated with *k*_{h} is

At larger *k*_{h}, on the other hand, Ro(*k*_{h}) > 1 and the dynamics will not be QG. If the dynamics at these *k*_{h} are stratified turbulence, the associated vertical wavenumber will be the buoyancy wavenumber *N*/*U* (Billant and Chomaz 2001; Waite and Bartello 2004; Lindborg 2006). Using the energy spectrum to define a *k*_{h}-dependent *U*, the buoyancy wavenumber at each *k*_{h} is

It is straightforward to show that *m*_{QG}(*k*_{h}) = *m*_{b}(*k*_{h}) when Ro(*k*_{h}) = 1. Assuming a mesoscale energy spectrum of the form *E*(*k*_{h}) = *cε*^{2/3}(*k*_{h}Δ*k*)^{−5/3}, where *ε* is the kinetic energy dissipation rate and *c* is an *O*(1) constant (e.g., Riley and Lindborg 2008), the expression for *m*_{b}(*k*_{h}) becomes

which is only weakly dependent on *k*_{h}.

The QG and buoyancy vertical wavenumbers are plotted in Fig. 4b for the highest-resolution simulation. They cross at *k*_{h} = 52, corresponding to a wavelength of 77 km, at which Ro(*k*_{h}) = 1. For *k*_{h} < 52, Ro(*k*_{h}) < 1 and the QG wavenumber is smaller than the buoyancy wavenumber. For these scales, *m*_{QG} < 26. For *k*_{h} > 52, the buoyancy wavenumber is smaller than the QG wavenumber. For these scales, *m*_{b} < 26 again (apart from the dissipation range at very large *k*_{h}, where *m*_{b} grows because of the drop-off in kinetic energy). As a result, to resolve the QG and stratified turbulence vertical scales associated with the full range of *k*_{h} outside the dissipation range, we need *m*_{T} > 26. Only the two highest-resolution simulations with *n*_{υ} = 128 and 256 satisfy this criterion, which are in fact the only two simulations that seem well converged. The two converged simulations in Fig. 3 have a mesoscale transition wavenumber *k*_{m} ≈ 25, at which Ro(*k*_{m}) ≈ 0.6.

The importance of (21) in underresolved cases is further investigated by considering additional simulations with *n*_{υ} = 16 and different values of *N* and *f*. Figure 5 shows energy spectra from four simulations in which *N* and *f* are doubled separately and together. The changes to *N* and *f* influence the level of large-scale kinetic energy, so spectra are compared after normalizing with the energy in *k*_{h} = 3 (Fig. 5b). When *N* is doubled relative to U16, the horizontal wavenumber given by (21) decreases by a factor of 2; the spectral transition obtained in Fig. 5 moves upscale consistently, and the resulting mesoscale spectrum is much more pronounced. Similarly, when *f* is doubled relative to U16, the transition wavenumber moves downscale. When *N* and *f* are both doubled, the transition also moves downscale slightly, though not as much as when *f* is doubled. These findings support the hypothesis that, in simulations that are very underresolved in the vertical, the transition to the mesoscale spectrum is determined by the vertical resolution through (21).

The strength of the kinetic energy cascade can be measured by the spectral energy flux

where is the nonlinear transfer spectrum [i.e., the advective term in the evolution equation for *E*(*k*_{h})] (e.g., Augier and Lindborg 2013; Deusebio et al. 2013; Peng et al. 2015). The value of Π(*k*_{h}) is the rate at which kinetic energy is transferred downscale through wavenumber *k*_{h}. Negative values denote upscale transfer. The flux is plotted in Fig. 6 (red curves) for *n*_{υ} = 16, 64, and 256. The shape of the flux spectrum is similar at all resolutions: there is upscale flux at the smallest two wavenumbers and a broad plateau of constant flux over 6 100. This constant flux range corresponds to a downscale cascade of kinetic energy through the mesoscale. The strength of the cascade is slightly lower in the lowest resolution case.

The shallowing of the mesoscale energy spectrum is frequently interpreted in terms of the different dynamics of balanced vortices and inertia–gravity waves (e.g., Bartello 1995; Hamilton et al. 2008; Deusebio et al. 2013). Figure 7 shows spectra of geostrophic and ageostrophic energy from simulations with the lowest, intermediate, and highest vertical resolution. These spectra are computed with the geostrophic/ageostrophic normal mode decomposition of the Boussinesq equations (Bartello 1995), which is a linear decomposition into geostrophic and wave modes. The geostrophic and ageostrophic energy sum to the total (kinetic plus available potential) energy. With low resolution (Fig. 7a), the geostrophic spectrum has a slope of −3 at large scales, below which it shallows significantly. For the ageostrophic spectrum, there is not a wide enough spectral range to measure the large-scale slope, but the mesoscale slope is very shallow at around −0.5. With higher resolution (Figs. 7b,c), the geostrophic spectrum again has a slope of around −3 but, unlike in the low-resolution case, the −3 range extends deeper into the mesoscale. The ageostrophic spectrum has a broad shallow slope at larger scales (−1.7 for *n*_{υ} = 64 and −1.5 for *n*_{υ} = 256, measured over 6 ≤ *k*_{h} ≤ 20) that also extends into the mesoscale beyond, below which it shallows significantly.

In more comprehensive atmospheric models with nonconstant *N*, the geostrophic/ageostrophic normal mode decomposition is more complicated, and a Helmholtz decomposition of the horizontal velocity is often applied instead (e.g., Hamilton et al. 2008; Skamarock and Klemp 2008; Waite and Snyder 2009, 2013). We have verified that the spectra of rotational and divergent horizontal kinetic energy look very similar to the spectra of geostrophic and ageostrophic energy in Fig. 7a: the rotational spectrum is approximately equal to the geostrophic spectrum, and the divergent spectrum is approximately equal to half the ageostrophic spectrum (the factor of 2 results from the fact that the ageostrophic energy has some potential energy, while the divergent energy is purely kinetic). In fact, apart from the very shallow range at the largest *k*_{h}, our high-resolution spectra are similar to those in Waite and Snyder (2009), who found the rotational and divergent energy spectral slopes to be around −3 and −5/3, respectively. These distinct spectral slopes are consistent with a downscale cascade of gravity waves driven by nonlinear interactions with the geostrophic modes (Bartello 1995), possibly by propagation through the vortical straining flow (Plougonven and Snyder 2005; Waite and Snyder 2009).

For a physical picture of the dependence on vertical resolution, Fig. 8 shows horizontal slices of the vertical vorticity and vertical velocity components at *t* = 24 days from the simulations with *n*_{υ} = 16 and 256. While both plots exhibit some noise, it appears much smaller in scale in the run with higher vertical resolution. The high-resolution case also has more filamentary and spiral structures in the vorticity field extending to smaller scales than in the case with lower *n*_{υ}. The larger-scale noise with *n*_{υ} = 16 is consistent with the inconsistent representation of QG and stratified turbulence motions at intermediate and small horizontal scales.

### b. With vertical mixing

In practice, most models are run with some vertical mixing, either through numerical viscosity and diffusion and/or ad hoc mixing schemes. To evaluate the effect of vertical mixing on the energy spectrum, we have repeated the above simulations with vertical hyperviscosity and a stability-dependent eddy viscosity scheme.

#### 1) Hyperviscosity

Figure 9 shows the kinetic energy spectra for the simulations with vertical hyperviscosity at different vertical resolutions. These spectra also exhibit significant dependence on *n*_{υ} for the two most underresolved cases. However, by contrast to the simulations with no mixing, the mesoscale spectra here are suppressed rather than amplified when *n*_{υ} is too low. In fact, the transition wavenumber seems to be decreasing toward its converged value rather than increasing (Table 1), and the mesoscale spectra with low vertical resolution are spuriously steep. In particular, the kinetic energy in *k*_{h} = 40 increases by a factor of 4 as *n*_{υ} increases from 16 to 256. The spectra are well converged for *n*_{υ} ≥ 64.

When *n*_{υ} = 16, the hyperviscosity case has much less mesoscale energy than the case with no vertical mixing (Fig. 10a). Relative to their higher-resolution counterparts, it is clear that the low-resolution mesoscale spectrum in the hyperviscosity run is artificially steepened, while the spectrum in the run without mixing is artificially amplified. Correspondingly, with low vertical resolution the spectral flux through the mesoscale is dramatically reduced compared to the case with no vertical mixing (Fig. 6a). By contrast, when *n*_{υ} = 64, the spectrum with hyperviscosity agrees well with the spectra with no vertical mixing out to around *k*_{h} = 40, beyond which the hyperviscosity spectrum is much steeper (Fig. 10b). When *n*_{υ} = 256, the spectrum with hyperviscosity has the same shape but slightly lower amplitude than the spectrum with no vertical mixing (Fig. 10c). The spectral flux exhibits a broad range of constant flux in the *n*_{υ} = 64 and 256 cases, but its amplitude is lower than the flux obtained with no mixing (Figs. 6b,c).

#### 2) Stability-dependent vertical mixing

Energy spectra for the simulations with the stability-dependent eddy viscosity are shown in Fig. 11. As found in the cases with no mixing, the shallow mesoscale spectrum is exaggerated when the vertical resolution is low, and the transition wavenumber increases with increasing *n*_{υ}. In fact, the spectrum obtained with *n*_{υ} = 16 is almost identical to the spectrum obtained with no vertical mixing (Fig. 10a). However, the resemblance is lost at higher resolution. While *k*_{m} increases as *n*_{υ} increases from 16 to 32, the shallow mesoscale spectrum is diminished as *n*_{υ} increases to 64, and it essentially disappears for *n*_{υ} = 128 and 256, which are almost indistinguishable from one another. Unlike the simulations with no mixing and hyperviscosity, the stability-dependent scheme gives no spectral transition for *n*_{υ} ≥ 128 (Figs. 10 and 11). The convergence of the spectrum for *n*_{υ} ≥ 128 is probably facilitated by the fixed mixing length ℓ = 30 m in the eddy viscosity.

The behavior of the stability-dependent mixing cases can be understood by considering how the values of Ri vary with vertical resolution. Figure 12 shows probability distributions of Ri from these simulations. The fraction of the domain with Ri < 0.25, in which the vertical mixing is enhanced, is extremely dependent on *n*_{υ}, because larger *n*_{υ} allows for the resolution of stronger shear and, hence, smaller values of Ri. When *n*_{υ} = 16, only 0.91% of the domain has Ri < 0.25; as a result, the vertical eddy viscosity in (12) is negligible almost everywhere and the mixing scheme has very little effect. As *n*_{υ} increases, this domain fraction also increases to 5.2%, 13%, 15%, and 14% for *n*_{υ} = 32, 64, 128, and 256. As a result, the vertical mixing is much more active in the high-resolution cases.

This mixing has a significant influence on the flux of kinetic energy through the mesoscale, as seen in Fig. 6. When the vertical resolution is low, the mixing is weak and there is a broad range of constant spectral flux, almost identical to the case without mixing (Fig. 6a). However, at higher vertical resolution, the strong mixing reduces the spectral flux over a wide range of *k*_{h}, thereby suppressing the cascade of energy to small horizontal scales (Figs. 6b,c). For *n*_{υ} = 64, this reduced flux is still larger than that obtained with hyperviscosity (Fig. 6b). However, for *n*_{υ} = 128 (not shown) and 256 (Fig. 6c), the flux is weaker than that obtained with hyperviscosity for *k*_{h } 10.

## 4. Discussion

High-resolution atmospheric simulations have dynamical structures with a wide range of vertical length scales, from deep QG vortices to thin shear layers. To accurately simulate these features, the vertical grid spacing must be small enough to resolve these vertical length scales. Since the vertical scale of structures generally decreases with decreasing horizontal scale, the smallest vertical length scales will be dependent on horizontal resolution. This paper has investigated the implications for the horizontal wavenumber energy spectrum of not resolving this full range of vertical scales. Simulations were analyzed of a forced baroclinically unstable jet with a wide range of Δ*z* from 2 km down to 100 m. With high vertical resolution that resolves the buoyancy scale, the energy spectra are converged and are insensitive to further increases in vertical resolution. However, at lower vertical resolutions, the spectrum is found to be very dependent on resolution. The nature of the vertical mixing scheme has a profound effect on this dependence.

With no vertical mixing or with weak, stability-dependent vertical mixing, the shallow mesoscale part of the energy spectrum is amplified when the vertical resolution is too low. The shallow part of the spectrum has a wider spectral range and a shallower slope than in the well-resolved cases. This behavior is observed here for the two lowest-resolution cases *n*_{υ} = 16 and 32, corresponding to Δ*z* = 2 and 1 km. While Δ*z* = 2 km is extremely coarse, 1-km-level spacing is not uncommon in the upper troposphere (e.g., Hamilton et al. 2008). These grids do not resolve the mesoscale buoyancy scale in the vertical direction and do not resolve the full range of QG vertical scales at larger horizontal scales with Ro < 1.

The exaggerated shallowing of the mesoscale spectrum in these cases appears to be due to the inconsistent grid aspect ratio of motions with small horizontal scales, as discussed by Lindzen and Fox-Rabinovitz (1989). As energy is transferred downscale in the horizontal, it is also transferred to small vertical scales. In the underresolved cases, there is a certain horizontal scale that corresponds to the smallest resolved vertical scale. As energy cascades further downscale, it is trapped at this finite vertical grid scale, which is artificially large for these small horizontal scales. These motions are forced by the grid to have an aspect ratio that is inconsistent with the underlying dynamics, and are therefore not represented correctly. There is no reason to expect this poorly resolved motion to be balanced, and unconstrained by potential vorticity conservation, it will be able to cascade efficiently downscale. Indeed, our underresolved simulations have a robust downscale energy cascade, which is driven in part by this spurious vertical-grid-scale motion.

On the other hand, at low resolution with vertical mixing (e.g., hyperviscosity), energy that makes it into the smallest resolved vertical scales is damped, not trapped. As a result, energy is dissipated at intermediate and small horizontal scales associated with vertical grid- and subgrid scales. In this case, the energy spectrum is unrealistically steep, and the downscale cascade of energy and transition to a shallow mesoscale spectrum is suppressed.

With high vertical resolution, the energy spectra obtained in simulations with no vertical mixing and with vertical hyperviscosity seem to converge. The converged spectrum exhibits the characteristic shallowing in the mesoscale with a transition scale of around 160 km, which is smaller than that in the exaggerated spectra with low resolution. Reasonable convergence is found in simulations with *n*_{υ} = 128 and 256, corresponding to Δ*z* ≈ 100 and 200 m. The converged spectra are similar in many ways to the Nastrom and Gage (1985) spectrum and spectra from other computational studies, but there are also some important differences. The transition scale of 160 km is in line with values reported elsewhere (e.g., Cho and Lindborg 2001); the transition is associated with the crossing of the steeper geostrophic spectrum and shallower ageostrophic spectrum; and both geostrophic and ageostrophic spectra shallow in the mesoscale (e.g., Lindborg 2007; Callies et al. 2014). On the other hand, the spectral slopes at small horizontal scales are much shallower than −5/3. The physical space plots suggest that these scales are dominated by small-scale stratified turbulence, more so than in idealized dry baroclinic wave simulations, in which inertia–gravity waves dominate at small scales (Waite and Snyder 2009, 2013). The present simulations support the hypothesis that stratified turbulence might partially explain the shallowing of the mesoscale spectrum, but more work is required to investigate the discrepancies between different studies, for example, by considering different initial jets and forcings. It is possible that greater horizontal resolution could yield a wider stratified turbulence inertial range with a slope closer to −5/3.

When a stability-dependent eddy viscosity is used for vertical mixing, the transition to the shallow mesoscale spectrum is suppressed at high vertical resolutions. In these cases, the fine vertical grids resolve significant regions of small Ri, leading to strong vertical mixing, which damps a wide range of horizontal scales. These results suggest that, at least for the purposes of reproducing the atmospheric kinetic energy spectrum, stability-dependent mixing schemes may not be worth the additional expense: at low resolution the scheme does very little, and at high resolution it is far too dissipative compared to hyperviscosity. It is possible that different results would be obtained with a different vertical mixing scheme (e.g., one with a more realistic stability function or a smaller mixing length ℓ at small Δ*z*; note, however, that even our highest vertical resolution had Δ*z* > ℓ). Recent work on large eddy simulation of stratified turbulence has shown that eddy viscosity models without any explicit Ri dependence work well as long as *L*_{b} is sufficiently resolved (Khani and Waite 2014, 2015).

These findings raise the troubling possibility that some model energy spectra may be artifacts of coarse vertical grids. Such spectra may look qualitatively correct, but may actually transition to a shallow mesoscale power law at larger scales that they should. Indeed, our spectra were only converged for Δ*z * 200 m. Such fine vertical grids are not commonly employed outside the boundary layer. These results indicate that computational studies of the energy spectrum must consider sensitivity to Δ*z* along with other parameters (as in e.g., Brune and Becker 2013). As argued by Lindzen and Fox-Rabinovitz (1989), increases in computational power should be invested toward increasing both horizontal and vertical resolution consistently. There may be little value in increasing the horizontal resolution if the newly resolved small horizontal scales are contaminated by an inconsistent grid aspect ratio.

It is interesting to note that, while the amplitude and slope of the mesoscale spectrum may be very sensitive to vertical resolution, other spectral quantites are more robust. In particular, the spectral flux and the relative amounts geostrophic and ageostrophic energy are not overly sensitive to resolution, at least when vertical mixing is not too strong. Even at low vertical resolution, simulations are able to reproduce a wide spectral range of constant energy flux, suggestive of a mesocale cascade. At all resolutions, the geostrophic energy spectrum is relatively steep and the ageostrophic spectrum is relatively shallow, and the transition is associated with the crossing of these spectra (however, the length scale of this transition is very dependent on Δ*z*). While this robustness is encouraging, it implies that these features should not be used as evidence of the correctness of a simulation. Just because a model is able to simulate a wide range of constant downscale energy flux does not mean that the cascade is being modeled correctly. Sensitivity to vertical resolution must always be checked directly.

## Acknowledgments

Model data used in this study are available upon request to the author. This paper benefited from comments by Sina Khani, Chris Snyder, and three anonymous reviewers. Karthik Velakur assisted with code optimization. Some computations were performed on facilities of the Shared Hierarchical Academic Research Computing Network (SHARCNET: http://www.sharcnet.ca) and Compute/Calcul Canada. Financial support from the Natural Sciences and Engineering Research Council of Canada is gratefully acknowledged.

## REFERENCES

*IUTAN Symposium on Turbulence in the Atmosphere and Oceans*, D. G. Dritschel, Ed., Vol. IUTAM Bookseries 28, Springer, 117–130.

^{−5/3}law inertial range in mesoscale two-dimensional turbulence