## Abstract

A new version of the dynamic Smagorinsky model is presented that applies for nonisotropic momentum diffusion in high-resolution atmospheric circulation models. While the horizontal mixing length is computed in accordance with scale invariance in the mesoscale regime of the horizontal energy cascade, the associated dynamic vertical mixing length (DVML) is based on a recently developed scale invariance criterion and represents an application of the scaling laws of stratified macroturbulence. The DVML is validated in high-resolution simulations with the Kühlungsborn mechanistic general circulation model, using triangular spectral truncation at wavenumber 330 and a vertical level spacing of about 200 m in the upper troposphere. For a proper choice of the test filter, the model simulates a realistic horizontal kinetic energy spectrum in the troposphere along with a realistic intensity of the Lorenz energy cycle. This result is obtained without any hyperdiffusion, and it depends only little on whether the vertical mixing length is prescribed or set to the DVML. The globally averaged Smagorinsky parameter is about *c*_{S} ≅ 0.53. The latitude–height cross sections show that *c*_{S} maximizes in regions of strong mesoscale kinetic energy.

## 1. Introduction

In the free atmosphere, the transfer of kinetic energy from large scales (baroclinic Rossby waves) to small-scale turbulence corresponds to a forward horizontal energy cascade that gives rise to a macroturbulent horizontal kinetic energy spectrum in the mesoscales (e.g., Nastrom and Gage 1985; Lindborg 2006; Hamilton et al. 2008; Brune and Becker 2013; Augier and Lindborg 2013), known as the Nastrom–Gage spectrum. General circulation models (GCMs) with high resolution are typically truncated at a horizontal resolution scale of about 10–100 km. Some explicit or implicit parameterization of horizontal momentum diffusion must be introduced to balance the forward energy cascade. Though this measure is often regarded mainly as part of the dynamical core of a model (i.e., as part of the applied algorithms to allow for numerically stable integration of the primitive equations), we consider horizontal diffusion primarily as an explicit parameterization of subgrid-scale processes (e.g., Meneveau and Katz 2000; Becker 2001; Lesieur et al. 2005).

The usual approach to balance the forward energy cascade in the mesoscales is to parameterize horizontal momentum diffusion by a linear biharmonic or higher-order hyperdiffusion (e.g., Collins et al. 2004; Stevens et al. 2013). Harmonic horizontal diffusion is commonly disregarded since it is not very scale selective in its linear form. Hyperdiffusion schemes, on the other hand, are sufficiently scale sensitive. They can be adjusted such that those scales that are not affected directly by the diffusion form a realistic horizontal kinetic energy spectrum (Takahashi et al. 2006; Siskind 2014). This is achieved by removing the kinetic energy that is injected into the spectrum at larger scales (e.g., baroclinic instability) below a transition scale that marks the so-called *effective resolution* (Skamarock 2004; Burgess et al. 2013). As stated by Skamarock (2004, p. 3027), “modes with wavenumbers higher than the effective resolution (modes with shorter wavelengths) are damped relative to their atmospheric counterparts and hence are dynamically suspect.” In this context we note that hyperdiffusion schemes are not based on a physically motivated turbulence model. Moreover, a viscous subrange forms that looks similar to the molecular viscosity subrange, albeit at much larger wavelengths, and often shows a buildup of energy at the smallest scales (e.g., Skamarock 2004; Hamilton et al. 2008). Hence, small-scale flow features like gravity waves cannot be studied below the effective resolutions scale.

The present study rests on the insight that, in a statistical sense, the mesoscales in the free troposphere represent a macroturbulent inertial range with a forward horizontal energy cascade. Therefore, we attempt to parameterize the subgrid scales in analogy to turbulence models used in large-eddy simulations (LES). To this end, two particular aspects must be taken into account. First, the resolution scale in common LES typically lies in the regime of isotropic Kolmogorov turbulence (i.e., the resolution scale is smaller than the Ozmidov scale, which is typically 10 m in the troposphere). On the other hand, the inertial range found in the atmospheric mesoscales is governed by the scaling laws of stratified turbulence (e.g., Lindborg 2006; Augier and Lindborg 2013), and the resolution scale of a high-resolution GCM (~50 km) lies in this regime. Therefore, a turbulence model for a GCM should take the scaling laws of stratified turbulence into account. Second, it is common practice in LES that the resolved −5/3 exponential spectral slope associated with the turbulent energy cascade continues down to the resolution scale. Hence, LES form a turbulent subrange without a pileup of the energy spectra within the resolved scales. Corresponding examples of LES can be found in simulations of turbulent channel flow (Park and Mahesh 2009), isotropic and stably stratified turbulence (Métais and Lesieur 1992; Khani and Waite 2015), turbulent boundary layer flow (Porté-Agel et al. 2000), and gravity wave breaking (Remmler et al. 2015). We shall follow this approach when formulating a turbulence model for high-resolutions GCMs. In other words, our goal is to simulate the Nastrom–Gage spectrum down to the resolution scale of the underlying GCM. Such a turbulence model would allow us to investigate small-scale features like gravity waves even close to the resolution scale. A few model studies exist where the simulated energy spectra show this behavior (Koshyk and Hamilton 2001; Kitsios et al. 2012).

Note that the parameterization of nonresolved scales by momentum diffusion and frictional heating (corresponding to the shear production of turbulent kinetic energy) represents a dissipative process. It parameterizes the cascade of kinetic energy toward unresolved scales, which is ultimately balanced by the dissipation of kinetic energy resulting from molecular viscosity (e.g., van Mieghem 1973). Technically speaking, because of this parameterized process, the strict definition of an inertial range being free from forcing and dissipation is not applicable here. However, we note that one has to distinguish between turbulent dissipation that represents the transfer of kinetic energy to smaller scales, and true viscous dissipation, which in the real atmosphere is ultimately due to molecular viscosity and violates scale invariance (Schaefer-Rolffs et al. 2015, hereafter SR15) In particular, in order to simulate the mesoscale macroturbulent inertial range in a forced-dissipative atmosphere, the subgrid-scale parameterization must always be dissipative.

We employ the classical Smagorinsky model (CSM) (Smagorinsky 1963, 1993) as a basis for the new turbulence model. The CSM can be considered as a generalization of the mixing-length concept based on the work of Prandtl (1942). When applied in an atmospheric GCM (Smagorinsky 1993; Becker and Burkhardt 2007) the strong anisotropy of horizontal and vertical scales results in a corresponding aspect ratio of the mixing lengths. The latter are prescribed in the CSM (e.g., as functions of height). Also, note that the diffusion is harmonic by definition, and that scale selectivity arises from the nonlinearity. The CSM represents the interaction between the resolved and unresolved scales in line with the hydrodynamic conservation laws and the second law of thermodynamics (Becker and Burkhardt 2007). However, the CSM is subject to two major problems: first, when using the CSM in a GCM with a nondissipative numerical core (e.g., without additional hyperdiffusion), a realistic Nastrom–Gage spectrum cannot be modeled by any choice of the mixing lengths (cf. Fig. 1); and second, the scale invariance constraint for the inertial subrange is violated, as was proven by Oberlack (2000) via Lie-group analysis. Note that scale invariance is equivalent to the presence of an inertial range (Kolmogorov 1941; Frisch 1995). Therefore, the Smagorinsky model must be modified in order to fulfill the scale invariance constraint.

A systematic mathematical approach to test whether a parameterization is scale invariant was recently given by SR15. Some implications of their scale-invariance criterion are noteworthy: first, the isotropic Smagorinsky model with a dynamically determined mixing length, the so-called dynamic Smagorinsky model (DSM; see Germano et al. 1991; Meneveau and Katz 2000), preserves scale invariance; and second, when the DSM is applied to the atmospheric mesoscales, scale invariance determines the horizontal mixing length, while the vertical mixing length depends on the horizontal mixing length according to the scaling laws of stratified macroturbulence. In this paper, we present a new anisotropic DSM that applies to the free atmosphere in high-resolution GCMs and complies with the concept of stratified macroturbulence. A major benchmark test of the new scheme is whether it allows us to simulate a realistic Nastrom–Gage spectrum down to the resolution scale.

The main idea of the DSM is to describe explicitly resolved structures with the same mixing-length ansatz that is supposed to represent the unresolved scales. Then, scale invariance can be used to determine the mixing length as function of the resolved flow. This approach has proven quite useful in technical fluid dynamics (Moin and Kim 1982) and recently in idealized simulations of stratified turbulence (Khani and Waite 2015). It has also been applied to LES of the atmospheric boundary layer for domain sizes of (1 km) (Porté-Agel et al. 2000; Kumar et al. 2006, 2010). However, it has been applied sparsely to the mesoscales in the free atmosphere. In a previous paper (Schaefer-Rolffs and Becker 2013, hereafter SRB13), we presented the first application of the DSM as a parameterization of horizontal momentum diffusion in a GCM with moderate horizontal resolution (spectral truncation at *n*_{T} = 120 corresponding to a horizontal resolution scale of ≃ 110 km) and driven by mechanistic differential heating. The application of the DSM resulted in a reasonable kinetic energy spectrum down to wavenumbers of about 90 (cf. Figs. 5 and 6 in SRB13). Yet, the transition to an inertial range with a −5/3 power-law spectrum in the mesoscales was not simulated, presumably because the resolution was still too coarse. Recently, we presented an improved DSM where the test filter used to determine the horizontal mixing length was separated from the resolution scale (Schaefer-Rolffs 2017). Simulations with a spectral truncation at *n*_{T} = 330 were promising regarding the Nastrom–Gage spectrum. However, a dynamic treatment of the vertical diffusion as part of the DSM has not been investigated yet. In combination with the horizontal DSM, such a turbulence model would represent a fully scale-invariant formulation of momentum diffusion in a highly resolved GCM. The present study is intended to propose and validate such a turbulence model.

The remainder of this paper is organized as follows. In section 2 we recapitulate the DSM for horizontal momentum diffusion. Section 3 introduces the calculation of the vertical mixing length according to the DSM for stratified turbulence. Numerical results are presented in section 4, which is followed by our conclusions (section 5).

## 2. Recapitulation of the dynamic Smagorinsky model

As mentioned in the introduction, the DSM for horizontal momentum diffusion is an extension of the corresponding mixing-length concept as described in Smagorinsky (1993) or Becker and Burkhardt (2007). The main difference in the DSM is that the horizontal mixing length *l*_{h} is no longer prescribed but is determined by the scale invariance constraint. The mixing length then becomes a function of the resolved flow according to

where is the resolution scale and *c*_{S} is the dynamically determined Smagorinsky parameter:

Here, a bar denotes quantities that are resolved in the model: is the resolved horizontal wind and is the trace-free horizontal shear tensor. Using indices 1 and 2 for zonal and meridional direction, the trace-free shear tensor can be written as

Here, and denote the resolved zonal and meridional wind, and denote the resolved horizontal divergence and relative horizontal vorticity, *ϕ* is latitude, and *a*_{e} is Earth's radius. The norm of the shear tensor is defined as || = ()^{1/2}.

We take advantage of the generalized DSM (gDSM; Schaefer-Rolffs 2017) to calculate *c*_{S}. The gDSM employs two test filter operations at two different filter scales , denoted by *tilde* and *roof* symbols, respectively. To find *c*_{S}, we define the stress tensors corresponding to the two filter scales as = − ( − ) and = −( − ). Since only the trace-free parts of the stress tensors are parameterized by the Smagorinsky model, we assume from now on and to represent the corresponding trace-free versions. Following SRB13 and Schaefer-Rolffs (2017), we compute *c*_{S} from the tensor equation:

by using a tensor-norm solution:

Here, the tensor is defined by the so-called Germano identity:

where the index 0 denote trace-free tensors. The tensor is given by

Then the turbulent diffusion coefficient is given by

and the turbulent stress tensor components become

Note that in SRB13 some approximations were used to simplify the and tensors [their Eqs. (26)]. For example, we neglected the cross-scale terms in our previous formulation, but found in the meantime that these terms are important. The formulation applied in this study includes the tensors and without any approximations. Consequently, the expressions used in Eqs. (6) and (7) are more accurate than the approximations used in SRB13.

## 3. The vertical mixing length according to the scale-invariance criterion

The scale invariance is usually not considered in the formulation of subgrid-scale parameterizations for GCMs or regional circulation models. SR15 derived a criterion that can be used to check analytically whether a parameterization is consistent with scale invariance. In the following we recapitulate this criterion and apply it to vertical momentum diffusion.

We define a scaling transformation of a quantity *a* using where *C*_{a} is the scaling parameter for *a*. SR15 showed that for horizontal momentum diffusion to be scale invariant, the horizontal mixing length *l*_{h} has to scale like the horizontal coordinate *x* (i.e., ). SR15 furthermore showed (see their appendix) that the scaling parameter for the vertical mixing length *l*_{z} differs from the scaling parameter *C*_{z} for the vertical coordinate *z*. Indeed, following SR15, scale invariance with regard to the horizontal scales is only fulfilled by the vertical momentum diffusion if the scaling parameters *C _{x}* and

*C*

_{z}for the horizontal and vertical coordinates obey a particular relation:

Hence, a functional dependence *C*_{z} = *f*(*C*_{z}) has to be sought for a unique solution of the problem. For this purpose we may exploit the scale-dependent aspect ratio for stratified turbulence (e.g., Lindborg 2006):

where *N* is the buoyancy frequency; *ε* is the dissipation rate; and *X* and *Z* are the horizontal and vertical length scales of the turbulence, respectively. Because *N* and are assumed to be constant, it follows *Z ***~ ***X *^{1/3}. There are two possibilities to combine Eqs. (10) and (11):

Equation (11) is valid for the vertical length scale (i.e.,

*C*_{z}=*C*_{x}/3). Then, Eq. (10) yields = 0. This corresponds to a constant (prescribed) VML, which is the usual approach for the free atmosphere in a GCM.- Equation (11) is valid for the vertical mixing length (i.e., = /3 =
*C*_{x}/3). In this case we get where*k*is a constant length scale. This parameter may be tuned such that*l*_{z}yields an averaged value that is typical for the prescribed asymptotic vertical mixing length in the free atmosphere as used in models.

At this point we cannot decide which result applies to the horizontal energy cascade in the mesoscales as both agree with the scale invariance criterion. A constant asymptotic vertical mixing length above the boundary layer is routinely used in our model. Results will be compared to those obtained in the second case where a dynamic vertical mixing length (DVML) according to Eq. (12) is applied. Note that we do not need to filter vertically in either case for the vertical diffusion to be with scale invariant. In particular, the traditional ansatz for the turbulent vertical diffusion coefficient in the free atmosphere,

is widely applicable, regardless of *l*_{z} being prescribed or given by Eq. (12).

## 4. Numerical results using a GCM with high resolution

### a. Model description

We apply a dry version of the Kühlungsborn Mechanistic general Circulation Model (KMCM). The KMCM is based on the spectral transform method and uses a terrain-following vertical hybrid coordinate according to Simmons and Burridge (1981). Except for the diffusion scheme, we use basically the same model setup that was used in Brune and Becker (2013) (i.e., T330 spectral resolution and 100 full model layers up to the stratopause, resulting in a level spacing of about 200 m in the upper troposphere and tropopause region). The horizontal grid consists of 1024 equidistant longitudes and 512 Gaussian latitudes. This ensures that all advection terms are transformed without spectral aliasing. The differential heating of the model is due to 1) temperature relaxation, 2) prescribed (and temporally constant) latent heating in the tropics, 3) self-induced latent heating in the extratropics, and 4) sensible heating resulting from the surface sensible heat flux. These methods are described in more detail in the appendix. The model includes an envelope orography (i.e., no orographic gravity waves are generated). This setup was also used in previous studies (e.g., Becker and Schmitz 2001; Körnich et al. 2006) and allows to simulate realistic dynamics for permanent January conditions. Such a concept of a mechanistic GCM works also well when extending the model vertically up the lower thermosphere (e.g., Becker 2009; Zülicke and Becker 2013; Becker and Brune 2014). More recently, the KMCM was equipped with explicit computations of radiative transfer and the tropospheric moisture budget [see Becker et al. (2015); Becker (2017), and references therein]. This new model version is currently used for studies of resolved gravity waves in the middle atmosphere (Becker and Vadas 2018). Some documentation of the zonal-mean climatology of the present version of the KMCM is given in the appendix. In the following we describe the details of our implementation of the gDSM.

The resolved scales subject to the test filter are expressed in terms of their total horizontal wavenumber. We specify the test-filter range by (90/120//120/150), where the first pair of wavenumbers indicates the wavenumbers where the filter sets in, while the second pair of numbers indicates where the filter turns off. The filter full width half maximum ranges from wavenumber 105 to 135. This formalism is described in more detail in SRB13 and in Schaefer-Rolffs (2017). Note that our test-filter regime lies in the regime of scales where we expect the kink in the Nastrom–Gage spectrum toward the −5/3 spectral lower law. This particular choice of the test filter will be further discussed later on. Our model runs were preceded by spinup integrations until an equilibrated kinetic energy spectrum was reached. The model was then integrated for another 18 days where every 45-min data were stored as snapshots. This time series forms the basis of the results presented below. We note that spectra of individual snapshot are always very variable, especially for the highest wavenumbers. This is actually caused by simulating some kind of backscattering of energy within the resolved scales, as discussed in Brune and Becker (2013), even though the averaged spectral energy flux is always downscale. Hence, all spectra shown hereafter are time averaged over one day.

The dynamically determined horizontal mixing length in the Smagorinsky scheme is given by [see Eq. (1)], where the Smagorinsky parameter is computed according to Eq. (5) and the resolution scale is given by (Orszag 1971)

Here, *n*_{T} = 330 is the truncation wavenumber, leading to .

The computation of horizontal diffusion tendencies is basically the same as in SRB13. The only exceptions are that replacements like → are necessary wherever tilde-filtered products emerge, and that the second term in the gradient of the Smagorinsky parameter,

is different from Eq. (29) in SRB13. In particular, the gradient of is computed according to^{1}

The implementation of the DVML is straightforward. In general, the vertical mixing length is given by

where *k*_{a} is the von Kármán constant and *z* is the actual height above ground. In the dynamic approach, the asymptotic vertical mixing length for the free atmosphere in Eq. (17), *l*_{asym} = 30 m, is replaced with *l*_{asym,dyn} = [cf. Eq. (12)]. The adjustable parameter is chosen to be *k* = 1.46 m such that the globally averaged DVML is about 36 m. Otherwise, the computation of vertical diffusion tendencies, including the local boundary layer scheme (Holtslag and Boville 1993), is the same as in SRB13 [see Becker (2003) for details of the discretization].

As in the model setup used in Brune and Becker (2013) and in SRB13, the present model also includes the full frictional heating owing to horizontal and vertical momentum diffusion according to the hydrodynamic energy conservation law. We also include a sponge layer in the stratopause region above 30 hPa, which is provided by a linear harmonic horizontal diffusion. Note that this particular sponge layer preserves angular momentum and total energy. Furthermore, we use the dynamic diffusion coefficients that are computed for momentum and also for sensible heat, assuming horizontal and vertical Prandtl numbers of 2 and 1, respectively. These Prandtl numbers were also used in the previous studies.

### b. Robustness of simulation results

Since scale-selective damping of the horizontal motion is required for numerically stable integration of a climate model, the new gDSM needs to be validated in this respect. As in the case of the CSM, model integrations are also stable using the new gDSM. As discussed in SRB13, numerical stability using the DSM mainly derives from the tensor norm solution of the Germano identity [see Eqs. (4) and (5)]. In particular, numerical stability requires that the diffusion coefficients are positive definite. This is also required by the second law since the macroturbulent diffusion scheme is an internal process of our model atmosphere (i.e., the associated frictional heating rates must be positive definite). This constraint would be violated if negative diffusion coefficients were allowed.

To demonstrate the numerical stability of the model integration we constructed a time series of daily mean energy spectra. Six spectra are shown in Fig. 2. The different colors correspond to different days of the spinup for the simulation with the gDSM including the DVML (red: day −14; orange: day −11; yellow: day −8; green: day −5; blue: day −2; violet: day 1 of the 18-day simulation used farther below). At synoptic scales (wavenumbers 15–60) the model is still adjusting and reaches somewhat higher energies for later days. On the other hand, the mesoscales show hardly any further adjustment from day −8 on. In particular, there is no accumulation of energy at the smallest scales. The model simulates the flattening of the spectrum in the mesoscales, as well as a slight drop-off of the energy close to the resolution scale. This proves that the gDSM simulates the scale-selective damping necessary for a roughly realistic energy spectrum quite well.

This conclusion is further emphasized by considering an arbitrary snapshot of the divergence of the horizontal wind (horizontal divergence) over the North Atlantic Ocean in the lower stratosphere at 100 hPa (Fig. 3). Since the large-scale horizontal wind is generally dominated by the geostrophic flow, a noisy behavior at the grid scale may be expected for the horizontal divergence, which is also a marker for mesoscale gravity waves. Because of the lack of an explicit moisture cycle in the present version of the KMCM, the highest mesoscale activity is induced by the breakdown of baroclinic Rossby waves and shows up in middle latitudes. The finescale structure that is visible in Fig. 3 reflects the instantaneous gravity wave activity. A similar wave activity at similar scales and with even higher amplitudes is also visible, for example, in the KANTO model [see studies of Watanabe et al. (2006) and Sato et al. (2012)]. Note that the KANTO model is a gravity-resolving spectral GCM with a T213 spectral resolution.

### c. Comparison of spectra

Figure 4 shows the simulated spectra at about 500 hPa (lower troposphere, blue), 250 hPa (upper troposphere, green), and 100 hPa (lower stratosphere, red), respectively, when using the gDSM with the DVML. The spectra in the lower and upper troposphere in the synoptic regime (wavenumber 15–100) show a very similar behavior, simulating a spectral slope that fits approximately to the expected −3 slope, and ending with a kink around wavenumber 120 (corresponding to a horizontal wavelength of about 330 km). In the lower stratosphere, a kink in the spectrum occurs at smaller scales around wavenumber 40 (similar to, e.g., Koshyk and Hamilton 2001; Bierdel et al. 2016). This shift in the kink marks the transition from the troposphere to the stratosphere as discussed by Burgess et al. (2013). In the mid- and upper stratosphere, our model does not describe the small-scale dynamics sufficiently well because of the increasing vertical level spacing above 75 hPa and the sponge layer above 30 hPa. Therefore, we do not discuss results for these altitudes.

Unlike the results presented in Burgess et al. (2013), where the effective resolution corresponds to a wavenumber of about 100, our simulations show an approximate −5/3 law in the mesoscales for wavenumbers from 130 to 300 in the upper troposphere. Only for the highest wavenumbers (larger than about 320) does the energy fall below the −5/3 law. In particular, the spectral energies do not increase with wavenumber close to the resolution scale. Note that such an unrealistic behavior is obtained when using the CSM (our Fig. 1).

In Fig. 5, we examine the components of the spectrum around 250 hPa. Different colors correspond to total (red), rotational (blue), and divergent (green) kinetic energy. Both components have a transition wavenumber similar to the total energy: the rotational part at wavenumber 140 from a −3 slope to a −5/3 slope, and the divergent part at wavenumber 40 from a −3 slope to a slope somewhat shallower than −5/3. While the rotational energy dominates in the synoptic scales, the divergent energy (green curves) dominates in the mesoscales, with an intersection of both curves around wavenumber 140, at a somewhat smaller wavenumber than obtained by Burgess et al. (2013). This may be attributed to the damping of the horizontal divergence in their version of the ECMWF model. Figures 4 and 5 also include the total spectral energies that are obtained when using a constant asymptotic vertical mixing length of 30 m (black dotted curves). The model produces basically the same results in this case. Interestingly, the minor differences manifest differently: at larger scales, the spectral energy energies obtained with constant vertical mixing length are somewhat lower than with the DVML, while in the mesoscales the ratio is reversed. These differences are presumably caused by a slightly stronger energy cascade when using the DVML.

Summarizing, both configurations, the application of the gDSM in combination with the DVML or with a constant vertical mixing length, allow for a reasonable simulation of the global energy spectrum at all heights. A reasonable behavior in the mesoscales is obtained without invoking artificial hyperdiffusion or numerical dissipation. A nearly realistic spectral slope extends close to the resolution scale. Hence, the new gDSM allows us to resolve the macroturbulence in the mesoscales of the free atmosphere with an effective numerical resolution close to the numerical resolution scale. This is analogous to LES of flows subject to isotropic turbulence as discussed in the introduction.

### d. Latitude–height cross sections

In the following we illustrate the performance of the gDSM with the DVML in terms of latitude–height cross sections of zonally and temporally averaged fields. Figure 6 shows the horizontal mixing length (colors) together with the total vertical mixing length (contours) according to Eq. (17). For the vertical mixing length we note the transition from the prescribed boundary layer mixing length given by *k*_{a}*z* to the asymptotic part in the free atmosphere given by Eq. (12). Since the DVML is proportional to the cube root of the horizontal mixing length, it is clear that the maxima of *l*_{h} and *l*_{z} are located in the same regions. Overall, both mixing lengths show a similar structure in the free atmosphere. Because of the cube root dependence of *l*_{z} on *l*_{h}, the relative variations of *l*_{z} are much smaller than the relative variations of *l*_{h}.

Figure 7 shows the simulated zonal wind (contours) and the total eddy kinetic energy (i.e., the kinetic energy contained in the deviation from the zonal-mean flow, shaded) for the simulation with the DVML. Both quantities are simulated quite reasonably. In particular, the eddy kinetic energy is comparable to other model results and ECMWF analyses for January (e.g., Roeckner et al. 1996). This quantity is dominated by planetary-scale and synoptic-scale Rossby waves.

Since the Smagorinsky parameter *c*_{S} is determined mainly from the scales between and [cf. Eq. (6)], we expect a correlation between *c*_{S} and the mesoscale kinetic energy. Figure 8 compares the eddy kinetic energy for wavenumbers *n* ≥ 90 only (color shading) with the Smagorinsky parameter *c*_{S} (contours). Indeed, there exists a clear correlation in the troposphere with maxima around 400 hPa in the extratropics. The tropospheric maximum is stronger in the winter hemisphere because of stronger baroclinic waves and hence a stronger energy cascade in the mesoscales. There are also maxima in the extratropical lower stratosphere from about 100 up to 50 hPa in summer and up to 20 hPa in winter. These maxima presumably reflect the dissipation of upward-propagating inertia gravity waves. Again, the wintertime maximum is stronger. There are two distinct maxima in the tropical stratosphere near the equator. This feature could be a consequence of inertial instability of equatorial waves.

The link between the mesoscale kinetic energy and *c*_{S} is also supported by Fig. 9, which displays the Smagorinsky parameter *c*_{S} (color shading) and the horizontal diffusion coefficient *K*_{h} according to Eq. (8) (contours). Both quantities show quite similar patterns. Hence, in contrast to the CSM, where the variability of the diffusion coefficient is determined by the norm of the shear tensor, the variability is strongly controlled by the variability of *c*_{S} in the gDSM. This underscores the central importance of scale variance.

The mass-weighted global-mean Smagorinsky parameter is *c*_{S} ≃ 0.538. This value exceeds results for isotropic turbulence by Schumann (1993), Vreman et al. (1997), or Khani and Waite (2015) by a factor of 2. The difference can probably be attributed to the different turbulent regimes that are parameterized by the DSM [i.e., isotropic Kolmogorov turbulence in these former studies vs stratified macroturbulence in our case or to the application of the tensor-norm approach in (5) instead of a least squares approach, cf. the discussion in Schaefer-Rolffs (2017)].

The physical consistence of the horizontal and vertical momentum diffusion computed with the gDSM is further demonstrated by diagnosing the intensity of the Lorenz energy cycle (cf. chapter V in Lorenz 1967). For this purpose, we compute the global-mean values of the vertically integrated frictional heating rates according to Becker (2003). The values for the two simulations (with prescribed or dynamical determined vertical mixing length) are given in Table 1. Noting that the dissipation is dominated by the surface friction, the value for the horizontal dissipation in relation to the total intensity of the energy cycle is quite significant because the horizontal diffusion contributes to the dissipation of kinetic energy mainly above the boundary layer. A value of around 0.34 W m^{−2} in relation to a total intensity of about 2.5 W m^{−2} confirms that horizontal diffusion should be regarded a physical parameterization. We also calculated the spectral contribution of horizontal and vertical dissipation in different wavenumber bins in the mesoscales. To avoid the influence from the atmospheric boundary layer or the sponge layer, we integrated only the dissipation in the upper troposphere from 500 to 100 hPa, where stratified turbulence is expected to be most pronounced. Table 2 summarizes the corresponding results, showing that even though the removal of mesoscale energy in the upper troposphere resulting from the DSM is present in the full mesoscale range, the smallest resolved scales (from wavenumber 270 on) provide half of the dissipation. Hence, the gDSM proves to be scale selective. Moreover, the increase of dissipation with increasing wavenumber is consistent with the theoretical calculation of the effective eddy viscosity from Kraichnan (1976) for an energy inertial range with a given resolution scale in order to keep the spectrum straight without a pileup of energy at the smallest resolved scales, mimicking the transfer of kinetic energy to the subgrid scales rather than the conversion into heat. Note that the spectral flux to the unresolved scales is represented by this dissipation.

## 5. Conclusions

The physical basis to apply some subgrid-scale turbulent diffusion in the free atmosphere of an atmospheric general circulation model (GCM) is still an unsolved problem. Usually, the boundary layer vertical diffusion extends over all heights using an asymptotic vertical mixing length. In addition, some scale-selective damping of the horizontal motion (horizontal diffusion) is necessary to balance the forward horizontal energy cascade. While horizontal diffusion is often regarded mainly as a numerical measure that is part of the dynamical core of a GCM (in practice the applied measures are usually hyperdiffusion or numerical filtering), we argue that horizontal diffusion is primarily a physical parameterization of unresolved scales. In particular, when the spatial resolution of a GCM covers a good part of the mesoscales, we can assume that the nonresolved scales in the free atmosphere are subject to a macroturbulent inertial range that is subject to the scaling laws of stratified turbulence (e.g., Lindborg 2006; Brune and Becker 2013; Augier and Lindborg 2013). This insight offers the perspective to invoke physically motivated methods to parameterize the unresolved mesoscales in terms on an anisotropic turbulence model.

In a previous work, we presented the dynamic Smagorinsky model (DSM) as a parameterization for horizontal momentum diffusion (Schaefer-Rolffs and Becker 2013). This method specifies the horizontal mixing length *l*_{h} in accordance with scale invariance. The vertical mixing length *l*_{z} was still treated conventionally (i.e., with a prescribed asymptotic value). This restriction was relaxed in the present study. More specifically, we used a recently developed scale invariance criterion (Schaefer-Rolffs et al. 2015) in order to find a formulation of the asymptotic part of *l*_{z} that complies with horizontal scale invariance in the free atmosphere. We found that the solution depends on whether the scaling relation of stratified turbulence is applied to the length scales or to the mixing lengths. The first possibility leads to a constant *l*_{z}, while a dynamic vertical mixing length (DVML) according to *l*_{z} ∝ is obtained in the second case.

We equipped the DSM with a more flexible test filter than was used in Schaefer-Rolffs and Becker (2013). According to Schaefer-Rolffs (2017) we refer to this new DSM as the generalized dynamic Smagorinsky model (gDSM). In particular, the small-scale limit of the test filter can now be separated from the resolution scale. Furthermore, we applied a much higher spatial resolution than in Schaefer-Rolffs and Becker (2013). To simulate a nearly realistic behavior of the horizontal kinetic energy spectrum we found it necessary that the test filter lies in the large-wavenumber regime of the resolved macroturbulent inertial range (i.e., at wavenumbers around 100) (corresponding to horizontal wavelengths of about 300 km), which are close to the kink in the energy spectrum. This finding can be interpreted as follows: the most reliable result for the dynamic computation of the Smagorinsky parameter from scale invariance [i.e., from the Germano identity, see Eq. (5)] seems to be obtained when the scales subject to the test filter are not affected directly by the diffusion. Hence, we assume that the test filter should lie at the low wavenumber end of the spectral range that is governed by the −5/3 power law. However, this interpretation needs to be confirmed in simulations with higher resolutions.

Test simulations with a dry version of a spectral GCM (equipped with mechanistic differential heating that allows us to simulate a realistic general circulation for permanent January conditions) show for the first time a nearly realistic simulation of the Nastrom–Gage spectrum without invoking hyperdiffusion or implicit numerical dissipation (Fig. 5). In particular, there is only a weak drop-off of the kinetic energy at the smallest resolved scales. Hence, the effective resolution of the model is very close to the numerical resolution scale. We also found that the dissipation due to horizontal momentum diffusion contributes quite significantly to the global mean dissipation. The comparison of simulations using a constant or dynamical vertical mixing length show only little difference. In particular, the scaling ratio *l*_{z} ∝ from stratified turbulence works well in the gDSM.

We obtained a global-mean value for Smagorinsky parameter of about 0.53. This value is about twice as strong as what is typically obtained when applying the DSM in LES of isotropic three-dimensional turbulence (e.g., Schumann 1993; Vreman et al. 1997; Khani and Waite 2015). Latitude–height cross sections confirmed that the spatial patterns of the Smagorinsky parameter *c*_{S} and the horizontal diffusion coefficient *K*_{h} are mainly determined by the mesoscale kinetic energy. The correspondence of *K*_{h} with *c*_{S}, rather than with the norm of the shear tensor as is the case for the CSM, indicates that introducing scale invariance to the Smagorinsky ansatz is essential. This conclusion is also evident when comparing the energy spectra obtained with the CSM (without any additional hyperdiffusion, see Fig. 1) to the results obtained with the gDSM (e.g., Figs. 2 and 5).

In the present study, we focused on scale invariance of momentum diffusion, and we applied the obtained horizontal and vertical diffusion coefficients to parameterize the anisotropic diffusion of sensible heat (using Prandtl numbers of 2 and 1, respectively). This approach still deviates from a strict realization of the concept of scale invariance. The reason is that stratified turbulence also implies a forward cascade of available potential energy that is comparable to the cascade of kinetic energy (Lindborg 2006; Augier and Lindborg 2013). Therefore, the thermal diffusion coefficients should be derived from the scale invariance constraint independently from the mechanical diffusion coefficients. A corresponding theoretical framework was proposed in Schaefer-Rolffs et al. (2015). In the future, we will include this concept in the gDSM. In this context we also note that the usual formulation of the turbulent vertical heat flux in terms of the vertical gradient of potential temperature has recently been questioned by Gassmann (2018) who suggested different formulations of the thermal diffusion for unstable and stable stratification such that the thermal dissipation is positive definite.

Another issue that needs to be addressed is the inclusion of a dependence of the diffusion coefficients on the Richardson number. Usually, the vertical diffusion coefficient in Eq. (13) is extended by multiplying with a so-called Richardson criterion (e.g., Holtslag and Boville 1993). This is the usual approach in a diagnostic turbulence model (i.e., a parameterization that does not include a prognostic equation of the turbulent kinetic energy, but instead presumes its quasi-stationary form) in order to take into account the effects shear production and buoyancy production of turbulent kinetic energy in vertical diffusion schemes. The particular Richardson criterion used in the KMCM above the Prandtl layer is *F*(Ri) = (1 − 18Ri)^{1/2} for Ri ≤ 0 and *F*(Ri) = 1/(1 + 9Ri) for Ri > 0. Becker [2009, see his Eq. (10)] found it necessary to extend also the horizontal diffusion from the CSM with such a functional dependence on the Richardson number in order to perform gravity wave resolving simulations in the mesosphere. In the future, we will use the criterion of Schaefer-Rolffs et al. (2015) to validate and, if required, reformulate the Richardson criterion in order to ensure that the scale invariance constraint is fulfilled when the diffusion coefficients are multiplied by a Richardson criterion.

The gDSM can furthermore be extended to include also the turbulent diffusion of specific humidity in the free atmosphere. This is motivated by the fact that in stratified turbulence, the variance of any tracer variance is subject to a forward cascade. In addition, Waite and Snyder (2013) showed that latent heating is a significant source of buoyancy production of mesoscale kinetic energy, corresponding to an injection of kinetic energy at the mesoscales. As mentioned in section 4a, a high-resolution moist version of the KMCM has recently been developed. Future investigations using this new model version will have to address the following questions: How is the parameterization of convection related to the concepts underlying the gDSM, and how will an explicit simulation of the moisture cycle modify the performance of the gDSM? Last but not least, a theoretical concept needs to be developed for matching the gDSM with the boundary layer scheme. In the present study we simply used the scale invariance constraint to specify the asymptotic vertical mixing length.

Summarizing, we presented the concept of a generalized DSM to parameterize the anisotropic turbulent diffusion in the free atmosphere of a mechanistic GCM that resolves part of the mesoscales. The model shows a nearly realistic global kinetic energy spectrum. The particular benefits of the new method are that the diffusion scheme is consistent with the hydrodynamic conservation laws and with the scale invariance constraint, and that the effective resolution of the model is close to the numerical resolution, allowing us to investigate dynamical processes even close to the resolution scale.

## Acknowledgments

We thank Rahel Knöpfel for many useful discussions. The valuable comments of two anonymous reviewers are grateful acknowledged. This paper is a contribution to the project T1 (Mesoscale energy cascades in the lower and middle atmosphere) of the Collaborative Research Centre TRR 181 “Energy Transfer in Atmosphere and Ocean” funded by the German Research Foundation.

### APPENDIX

#### Differential Heating and Zonal-Mean Climatology

The differential heating applied in the present version of the KMCM is comparable to the Held–Suarez benchmark test (Held and Suarez 1994) but is significantly more realistic because of the generation of quasi-stationary Rossby waves by orography and latent heating. The differential heating (or external diabatic heating) is given by simple representation of radiative and latent heating, plus the vertical diffusion of heat. The vertical integral of the latter is given by the surface sensible heat flux. The sum of radiative and latent heating is written as

where *T* is temperature, *T*_{e} is a zonally symmetric equilibrium temperature, *ω* is the pressure velocity, and *h* is the Heaviside step function. The relaxation time *τ* increases from 16 days in the troposphere to 40 days at 100 hPa and decrease to 7 days in the upper stratosphere and mesosphere (Dunkerton 1991). We furthermore apply = 40 hPa day^{−1}. Moist convective heating in the deep tropics is represented by an additional prescribed cumulus heating *Q*_{c} (Hou 1993) and latent heating in the middle latitudes is accounted for by self-induced condensational heating (Mak 1994) with a prescribed heating function *Q*_{m}. The patterns of *T*_{e}, *Q*_{c}, and *Q*_{m} correspond to January conditions as presented in Fig. A1. Empirical formulas for *T*_{e}, *Q*_{c}, and *Q*_{m} are given in Becker and Schmitz (2001, see their appendix C).

Land–sea contrasts are incorporated by including the heating functions in the definition for the surface temperature:

The factor 0.4 is chosen to obtain a realistic surface temperature in connection with *Q*_{c} and *Q*_{m}. Figure A2 shows the model orography and surface temperature according to Eq. (A2). Since the orography and the surface sensible heat flux are included, the mechanistic model yields semirealistic representation of synoptic and planetary Rossby waves for the Northern Hemispheric winter season (Körnich et al. 2006). In particular, Fig. A3 shows that the zonal-mean zonal wind and temperature, as well as the meridional eddy momentum and heat fluxes, are simulated realistically (e.g., cf. plots in Randel 1992; Roeckner et al. 1996).

## REFERENCES

*Turbulence: The Legacy of A. N. Kolmogorov*. Cambridge University Press, 296 pp.

*Large-Eddy Simulations of Turbulence.*Cambridge University Press, 232 pp., https://doi.org/10.1017/CBO9780511755507.

*Symmetrie, Invarianz und Selbstähnlichkeit in der Turbulenz (Symmetry, Invariance, and Self-Similarity in Turbulence; Habilitation)*. RWTH Aachen, 187 pp.

*Führer durch die Strömungslehre (Essentials of Fluid Mechanics)*. Fried. Vieweg & Sohn, 105–108 pp.

*Lecture Series: Introduction to the Modeling of Turbulence*, W. Kollmann, Ed., Von Karman Institute,

*Large Eddy Simulation of Complex Engineering and Geophysical Flows*, B. Galperin and St. A. Orszag, Eds., Cambridge University Press, 3–36, doi:.

*Atmospheric Energetics.*Clarendon Press, 306 pp.

## Footnotes

© 2018 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).