High-resolution large-eddy simulations of stably stratified flows : application to subkilometer-scale turbulence in the upper troposphere – lower stratosphere

Large-eddy simulations of stably stratified flows are carried out and analyzed using the mesoscale atmospheric model Méso-NH for applications to kilometerand subkilometer-scale turbulence in the in the upper troposphere–lower stratosphere. Different levels of turbulence are generated using a large-scale stochastic forcing technique that was especially devised to treat atmospheric stratified flows. The study focuses on the analysis of turbulence statistics, including mean quantities and energy spectra, as well as on a detailed description of flow topology. The impact of resolution is also discussed by decreasing the grid spacing to 2 m and increasing the number of grid points to 8× 109. Because of atmospheric stratification, turbulence is substantially anisotropic, and large elongated structures form in the horizontal directions, in accordance with theoretical analysis and spectral, direct numerical simulations of stably stratified flows. It is also found that the inertial range of horizontal kinetic energy spectrum, generally observed at scales larger than a few kilometers, is prolonged into the subkilometric range, down to the Ozmidov scales that obey isotropic Kolmogorov turbulence. This study shows the capability of atmospheric models like Méso-NH to represent the turbulence at subkilometer scales.


Introduction
The environmental impact of aviation is a problem of increasing concern among scientists and policymakers as commercial air traffic continues to grow (IPCC, 1999;Lee et al., 2009).Because of long residence time, low-background concentrations and large radiative sensitivity at cruise altitude, aircraft emissions can influence the chemical and physical state of the atmosphere and affect the radiative budget of Earth.According to recent evaluations, the induced cloudiness produced by aircraft emissions in the form of contrails and induced cirrus is among the most uncertain contributors to the global radiative forcing (Sausen et al., 2005;Lee et al., 2009).One important reason for this uncertainty is that, although aircraft emissions are (generally) known at the nozzle exit, their impact depends on their interactions with the atmosphere and should be evaluated at the scales resolved by general circulation models or chemistry transport models, i.e., 100 km or more for grids that cover the entire globe.Predicting the evolution of emissions during all their lifetime is then a crucial step to evaluate their environmental impact and to produce parameterizations for large-scale models.This requires the use of atmospheric models to represent the smallscale motions and the physicochemical processes occurring in the free atmosphere.In particular, the dispersion of aircraft plumes and the mixing of exhausts with ambient air depend on the properties of atmospheric turbulence in the upper troposphere-lower stratosphere (UTLS) where most of flight time is spent during cruise (Gerz et al., 1998).
Published by Copernicus Publications on behalf of the European Geosciences Union.
Reproducing the effects on aircraft plumes of atmospheric turbulence in the UTLS using computational tools remains a challenging problem (Sharman et al., 2012).The usual approach employed in the literature is to represent aircraft plumes by Gaussian plumes and to initialize turbulence using idealized flow fields satisfying the given model spectra (Dürbeck andGerz, 1995, 1996;Gerz et al., 1998).This procedure imposes a limitation on the duration of the simulation since turbulence is not sustained and eventually decays.However, the vast majority of direct numerical simulations (DNS) and large-eddy simulations (LES) of forced stratified turbulence reported in the literature are based on nondimensional spectral formulations of Navier-Stokes equations under the Boussinesq approximation.While this allows for a fine description of turbulent structures and statistics against which mesoscale models can be validated, the inclusion of atmospheric processes like microphysics or radiation can be more problematic in these formulations.By contrast, in the present work, the properties of turbulence in the UTLS at the subkilometer scale are studied using a mesoscale atmospheric model formulated in physical space and coupled to a stochastic forcing technique that allows the nonlinear cascade to develop fluctuations at smaller scales.The choices of the model are motivated by the necessity of dealing with a model that contains all the physics (microphysics, radiative transfer, etc.) needed for studying environmental problems that are relevant to aviation and that will be addressed in follow-up studies.
The statistical properties of turbulence in the UTLS were first studied in the seminal paper by Nastrom and Gage (1985), who analyzed observational data from instrumented commercial aircraft during the GASP (Global Atmospheric Sampling Program) campaign.They observed the kinetic energy spectra vary with wave number k as k −5/3 for scales of a few kilometers to a few hundred kilometers, and as k −3 for larger scales.Similar results were obtained by Lindborg (1999), who analyzed data from the MOZAIC (Measurements of Ozone Water Vapor by Airbus Aircraft) campaign.Ideal turbulent, stratified flows, which have characteristics similar to those encountered in the UTLS, have been reproduced since then by means of DNS at different Reynolds numbers with or without the inclusion of hyperviscosity (Métais and Lesieur, 1992;Kaltenbach et al., 1994;Riley and Lelong, 2000;Riley and de Bruyn Kops, 2003;Waite and Bartello, 2004;Lindborg, 2006;Brethouwer et al., 2007;Riley and Lindborg, 2008).The picture emerging from these studies is that an inertial range exists in the horizontal energy spectrum, which is an indication of a downscale transfer of energy from large to small eddies as in classical Kolmogorov turbulence.However, because of atmospheric stratification, turbulence differs from isotropic turbulence, which reflects for example in the presence of large, elongated horizontal flow structures.Indeed, the competition between turbulence and stable stratification that tends to compress eddies vertically is the clue that makes the sim-ulation of strongly stratified flows challenging as turbulence can be three-dimensional yet substantially anistropic (Waite and Bartello, 2004;Lindborg, 2006;Brethouwer et al., 2007;Waite, 2011;Kimura and Herring, 2012).The intermediate range of scales between a few kilometers down to a few tens of meters of the Ozmidov scale where turbulence can be considered isotropic, is more difficult to observe and to predict numerically because of the high resolutions required to capture the small-scale dynamics.While some recent measurements (Wroblewski et al., 2010) seem to support a prolongation of the k −5/3 slope, no definitive answer has been found so far to elucidate the characteristics of this range of the spectrum.This is particularly important for the dispersion of aircraft emissions because the transition from the aircraft wakecontrolled to the atmosphere-controlled dynamics occurs exactly at these scales (Dürbeck and Gerz, 1996;Paugam et al., 2010).
The objective of the present study is to analyze the accuracy of the proposed method and LES models based on explicit transport of turbulent kinetic energy to reproduce the properties of turbulence in the UTLS described above.To that end, LES of stably stratified flows are performed in an idealized atmosphere using computational grids with up to 2 m resolution and up to 2048 3 grid points, one of the largest LES attempted in this area.The paper is organized as follows.Section 2 describes the governing equations of the model, with emphasis on the forcing methodology used to force turbulence at larger scales.The results of the simulations are discussed in Sect. 3 including the flow statistics, flow topology and energy spectra for different levels of atmospheric turbulence, and some elements of comparison with observational analysis.Conclusions are drawn in Sect. 4.

Governing equations and numerical model
The simulations were carried out using Méso-NH (mesoscale non-hydrostatic), the atmospheric research model developed by Météo France and the Laboratoire d'Aérologie.The model is built around a dynamical core that is capable of simulating atmospheric motions ranging from the meso-alpha scales down to the micrometer scales; an ensemble of packages treating different physical processes in the atmosphere; a flexible file manager, and an ensemble of preprocessing tools to set up the initial conditions, either idealized or interpolated from meteorological analysis or forecasts (for details, see the online scientific documentation: http://mesonh.aero.obs-mip.fr/mesonh).Briefly, the model solves Navier-Stokes equations in the anelastic approximation (Lipps and Hemler, 1982).The "reference" thermodynamic state is a function of the altitude (the vertical coordinate z).The basic prognostic variables are momentum u = [u, v, w] and potential temperature θ = T / where T denotes temperature and ≡ (p 00 /p) R d /C p d the Exner function where p and p 00 are the local and a reference ground level pressure, while R d and C p d are the gas constant and specific heat of dry air, respectively (Lafore et al., 1998).Additional prognostic variables include mixing ratios of chemical species and water substances limited to water vapor mixing ratio r v in this study.The virtual potential temperature is then defined as , where R v is the gas constant of water vapor.Coriolis forces are not important at the scales of interest for this study.Phase changes and microphysical processes are not activated.In the LES approach, each variable φ is decomposed into a resolved or filtered part φ and an unresolved or subgrid-scale part φ with φ ≡ φ + φ .This procedure can be obtained by a convolution integral of the variable with filter functions that depend on filter widths.In the present LES formulation, the filter widths correspond to the grid spacings x, y, and z, so that filtering is equivalent to grid averaging.With these hypotheses and neglecting third-order turbulent fluxes, the governing equations read (Lafore et al., 1998;Cuxart et al., 2000): where the subscript 0 denotes the reference state; g is the acceleration due to gravity, is the (filtered) modified pressure, and the specific heat of moist air.The term f denotes a generic body force, which corresponds to the stochastic forcing in this study and is discussed in the next section.All Reynolds stresses u ⊗ u , subgrid-scale fluxes u θ , u r v , and subgridscale correlations r v θ , θ 2 , are modeled using a closure based on mixing length L = ( x y z) 1/3 (Deardorff, 1972;Redelsperger and Sommeria, 1981) and a prognostic equation for the turbulent kinetic energy, which represents the isotropic part of Reynolds stress tensor, as explained by Cuxart et al. (2000): where C 2 m = 0.2 and C = 0.7 are model constants while the last term in the rhs (right-hand side) is the turbulent dissipation rate of kinetic energy, The subgrid-scale correlations involving potential temperature are also reported for the sake of completeness: where C s = 4 and C θ = 1.2 are model constants, while is a diagonal matrix defined as in Redelsperger and Sommeria (1981) and Cuxart et al. (2000).Finally, it is worth noticing that the present formulation of the model does not account for the contribution of molecular viscosity to the shear stress tensor.This is justified for atmospheric high Reynolds numbers if the cutoff length (proportional to the grid size) is much larger than the Kolmogorov dissipative scale η ∼ ν 3 / k 1/4 .For air molecular viscosity ν ∼ 10 −5 m 2 s −1 and turbulence dissipation rates at the tropopause k ∼ 10 −6 − 10 −5 m 2 s −3 (Lindborg, 1999), the Kolmogorov scale η ∼ 10 −3 m, which is much smaller than the grid sizes employed in this study.

Turbulence forcing
In order to obtain a statistically stationary velocity field, a low wave number body force is applied to the momentum equations using the method originally developed by Eswaran and Pope (1988) in spectral space and reformulated by Paoli and Shariff (2009) in physical space.In this method, all modes with wave numbers within a sphere of radius k f are forced using stochastic processes that mimic the turbulence production at scales larger than the computational domain and allow the nonlinear cascade to determine fluctuations at smaller scales.Denoting by R {•} and I {•} the real and imaginary part of a complex number, the body force is represented as a finite Fourier series, where i = √ −1, k is the three-dimensional wave number, and f(k, t) is the amplitude of the complex Fourier mode, which is obtained from the divergence-free projection, where τ f and σ f are, respectively, the autocorrelation time and the standard deviation of the process, while N | t+dt t and M| t+dt t are vector-valued normally distributed random processes with zero mean and unit variance over the time interval [t, t + dt].The mode corresponding to k = 0 is not forced, so from Eq. ( 10) the net force acting on the flow is zero and no mean motion is created.Furthermore, in the present implementation of the method, both vortical and divergent modes are excited, and the vertical component of the forcing is set to zero.This implies that vertical motions can be attributed to the nonlinear interactions between horizontal and vertical momentum rather than to the direct effect of forcing.As discussed by Eswaran and Pope (1988), in the isotropic turbulence case the kinetic energy dissipation rate resulting from the forcing scales as the product σ 2 f τ f .This relation has been verified for anisotropic turbulence (not shown).Hence, in this study, the sensitivity of turbulence to the forcing scheme has been analyzed by varying σ f and keeping τ f fixed.

Numerics
For the momentum components, the spatial discretization relies on a fourth-order finite differences staggered scheme, while time discretization is based on a leapfrog scheme, except for the stochastic body force that is advanced in time using a first-order Euler scheme.The piecewise parabolic method (PPM) (Colella and Woodward, 1984) advection scheme is used for all scalar variables.Pressure is obtained by solving a Poisson equation that enforces the incompressibility condition Eq. ( 1) by means of a FFT algorithm, which was specially devised for treating domain decomposition in massively parallel architectures.The code uses MPI as the basic network communication protocol and has proven good scalability on up to 100, 000 parallel cores during various benchmark tests of extreme-scale computing in Europe and in the US (Pantillon et al., 2011).For this study, the simulations were run on 4096 cores of a BullX supercomputer for an overall cost of 7 million CPU hours.

Computational details
The computational domain is a cubic box with side lengths, L = 4 km, representing a portion of an idealized atmosphere (see Fig. 1).The reference altitude is 11 km, which is typical of long-haul flights operated by commercial airlines.Values of temperature, pressure and density at this altitude are set to T ∞ = 218 K, ρ ∞ = 0.388 kg m −3 , and p ∞ = 24 286 Pa, respectively, and the corresponding potential temperature is θ ∞ = 326.6K.The background atmospheric fields ρ 0 (z) and θ 0 (z) are obtained from the hydrostatic and thermodynamic relations with the Brunt Väisälä frequency set to The background vapor-mixing ratio is obtained from climatology of a midlatitude summer standard atmosphere (Clatchey et al., 1972) with r v ∞ = 0.001.The flow is initially at rest, and the spectral shell of the forcing (highest forced wave number) is where is the fundamental wave number.Periodic boundary conditions are imposed in the horizontal directions x and y, whereas open boundary conditions are used in the vertical direction z in combination with two buffer zones of size δ = 200 m at the top and bottom of the computational domain where the velocity is smoothed down to zero.This boundary condition indirectly leads to the periodicity of velocity and temperature fluctuations in the z direction.However, data in the buffer zones are excluded from post-processing analysis.

Flow statistics
In order to characterize the flow statistics, the volume average φ of a generic quantity φ is defined as where L x = L y = L, and L z = L − 2 δ.Note that by definition of grid-averaging it follows that φ ≡ φ .The horizontal average φ h is defined in a similarly way as In the following analysis, the explicit dependence of φ and φ h on t and z will be omitted for notational ease.

Velocity and kinetic energy
Setting φ = u in Eq. ( 16), observing that the flow is initially at rest and the forcing produces zero net acceleration of the flow, yields where the last two terms of Eq. ( 19) represent, respectively, the resolved and subgrid-scale contributions to the secondorder moment of the velocity field.Setting φ = e k and using Eq. ( 5) yields the volume average of the subgrid-scale turbulent kinetic energy, while setting φ = E k ≡ u • u 2 and using Eqs.( 19)-( 20) yields which represents the total kinetic energy of the flow.In strongly stratified flows, where large anisotropy exists between horizontal and vertical directions, the resolved velocity u is further decomposed using Eq. ( 17) into a horizontal mean: and a perturbation with respect to this mean: with u h and u t satisfying the identity u t ≡ u − u h = u − u = 0. Using u t instead of u allows for a more natural interpretation of turbulence statistics that pertain to scales smaller than those that are directly forced.As observed by Chung and Matheou (2012), this decomposition discounts the energy of the larger horizontal modes from the total kinetic energy budget, and is particularly useful in the present simulations where the forcing produces local distortions of the flow that can be assimilated to nonuniform shear in the vertical direction.These motions correspond to the vertically sheared horizontal flow (VSHF) modes that were first identified by Smith and Waleffe (2002).
in Eq. ( 16) yields which represents the mean kinetic energy associated to u t .
It is interesting to observe that applying the volume-average operator to u t • u t and using Eq. ( 23) yields which shows that the second-order moments of u and u t coincide only if the horizontal mean is zero at each vertical level -a condition that is strictly verified for an infinite domain in the horizontal direction.Substituting Eq. ( 25) into Eq.( 21) and using Eq. ( 24), the mean kinetic energy E k can be recast as where the first two terms in the rhs represent the (resolved and subgrid-scale) kinetic energy of "genuine" fine-scale turbulence, while the last term represents the energy of shearlike motions (Smith and Waleffe, 2002;Lindborg, 2006).Introducing the variances of the velocity components of u t , the resolved turbulent kinetic energy can be recast as while the root-mean square of horizontal turbulent fluctuations is defined by which corresponds to the large-scale (horizontal) vortices (Waite, 2011).Finally, the horizontal Froude number is defined as in Brethouwer et al. (2007):

Potential temperature and potential energy
As for Eq. ( 22), potential temperature can be split into a horizontal mean, where the last equality is due to the smooth gradient of background temperature, and a perturbation with respect to this mean, Following the same procedure used for kinetic energy, the variance of potential temperature can be split as where again the first two terms in the rhs represent the contributions due to the resolved and subgrid-scale fluctuations, respectively, while the last term represents the contribution of the mean gradient of background temperature.The potential energy E p is related to the potential temperature variance through where θ ref ≡ θ 0 has been used as a reference temperature (Lindborg, 2006).It is known (see for example Kurien et al. (2006) that the sum of total kinetic and potential energy is conserved in the Boussinesq approximation so that E p is in general a more useful quantity for the analysis of flow energetics.Using Eqs. ( 31) and ( 32) and integrating over the domain yields with the usual notations: where the last two terms represent the potential energy of resolved and subgrid-scale motions, respectively.

Evolution of turbulent statistics
Three levels of turbulence, denoted by labels "S", "M" and "W" for strong, moderate and weak turbulence, were analyzed by varying the forcing parameter σ f in Eqs. ( 12) and ( 13).The timescale of forcing is τ f = 33.6 s.For numerical integration of Uhlenbeck-Ornstein process τ f has to be larger than the time step t of the simulation to avoid excessive random noise, and smaller than the largest timescale of the flow that can be estimated here as τ large = U 2 / k to avoid large-scale drift of the flow.In all cases considered here t ranges between 0.3 and 1 s while τ large is on the order of 10 4 s so that both conditions are satisfied.The effects of resolution on turbulence characteristics were also analyzed by varying from 10 to 2 m, with the corresponding number of grid points varying from 400 3 to 2048 3 .The nomenclature used for the different runs and values of the main statistical properties are documented in Table 1. Figure 2 shows the mechanism of turbulence generation through the stochastic forcing method for all the simulated cases.The evolution of the resolved turbulent kinetic energy E kt is characterized by an initial rapid increase of 0 < t < t 1 with t 1 on the order of 2 to 3 hours.This is followed by a transient phase, t 1 < t < t 2 with t 2 generally lesser than 5 hours, until the energy reaches statistically stationary conditions in the sense that it does not increase or decrease monotonically but oscillates around a steady value.This value increases with the forcing level but is roughly independent on the resolution, which confirms that the turbulence model has no or limited impact on the mean resolved turbulence.The evolution of turbulent potential energy E pt follows trends similar to the those of kinetic energy with the ratio E pt / E kt varying between 0.12 and 0.18, in the same range observed in previous DNS, e.g., Lindborg (2006).The ratio R = e k /( E kt + e k ) of the subgrid-scale to the total turbulent kinetic energy is shown in Fig. 3.As expected, this ratio increases when decreasing the resolution and the turbulence forcing, which reflects the fact that the contribution of the subgrid-scale stress to the total Reynolds stress increases.However, in all cases R < 0.1, which meets the classical requirement proposed by Pope (2000) for "well resolved" LES that the resolved turbulent kinetic energy should be at least 80% of the total energy.The evolution of the mean dissipation rate is shown in Fig. 4. The initial condition for e k corresponds to the threshold value that is initially assigned to the subgrid-scale energy to trigger grid-scale perturbations and initiate the process of energy transfer from the large, forced scales to the small, dissipative scales.Once the transfer is fully sustained by the forcing, the turbulence model should be able to create its proper dynamics and e k should increase, indicating that energy is correctly dissipated at the grid-scale level.This condition is verified for all cases except at the lowest resolution for low and moderate turbulence (run W10 where k k (0) during all the simulation time), which indicates a slightly inaccurate representation of turbulence at the smallest resolved scales  for these cases.This point is further discussed in the next sections.It can be observed that k starts to increases at t t 1 , corresponding to the effective activation of the turbulence model, and then attains a statistically steady value.Figure 5 shows the rms (root mean square) of turbulent velocities for the selected case M04.The figure provides a first indication of the anisotropic character of stratified turbulence, with horizontal fluctuations σ u and σ v much larger than vertical fluctuations σ w (the same behavior is observed for all other cases).The anisotropy ratio is σ u /σ w 0.3-0.4depending on the cases and is in the range of available observational data (Nastrom and Gage, 1985).

Turbulent structures
The turbulent structures of the flow are visualized in Fig. 6 by means of snapshots of potential temperature fluctuations It is also interesting to observe that, as forcing is increased, the initially thin and patchy layers become ticker and more chaotic, as the consequence of the increase in potential of turbulent eddies to overturn against background stratification.These flow characteristics can be quantified by monitoring the evolution of the buoyancy scale, that is reported in Fig. 7 for all the simulated cases (note the factor 2π accounts for the transformation from wave numbers to wavelengths, see Waite and Bartello, 2004;Waite, 2011).The buoyancy scales characterizes the thickness of the shear layers in stratified turbulence (Waite, 2011) and is also associated with the zigzag instability of columnar vortices (Billant and Chomaz, 2000) and overturning of internal gravity waves (Waite and Bartello, 2006).As shown in Fig. 7, L b varies between 250 and 300 m depending on the cases, which is well above the grid sizes used for the simulations.This guarantees that the largest vertical structures of the flow are resolved.The ratio L b /L z < 0.1, which also guarantees that the flow is not vertically confined.In addi-tion to the buoyancy scale, in LES of stratified turbulence it would be beneficial to resolve the Ozmidov scale, where the horizontal average S 2 h is obtained by setting φ = S 2 in Eq. ( 17).The Richardson number is often associated to the stability of stratified shear flows: for example, it is known from linear stability analysis (Miles, 1961;Howard, 1961) that an initially laminar flow becomes unstable (eventually evolving into a turbulent flow) below the critical value Ri c = 0.25; furthermore, values of Ri < 1 have been consistently measured in experimental observations and numerical simulations (Riley and de Bruyn Kops, 2003;Brethouwer et al., 2007).This picture is consistent with the vertical profiles and the probability density function (PDF) of Ri that are reported in Figs. 8 and 9, respectively.As a general trend, the Richardson number decreases when increasing turbulence forcing as this can produce locally stronger shear (Riley and de Bruyn Kops, 2003; Kimura and Herring, 2012).In the high-resolution cases, = 2 m, Ri fluctuates around 0.1 and 0.15 for the weak and moderate forcing, respectively.However, in the low-resolution cases, = 10 m, these values fluctuate around Ri > 1 for the weak forcing, Ri 1 for the moderate forcing, and Ri 0.5, well above 1, for the strongest turbulence.This again suggests that with a resolution of = 10 m the strongest turbulence is well resolved, the moderate turbulence is barely resolved and the weak turbulence is not resolved.To further verify this point, we compared the Ri profiles obtained from the 10 m resolution LES data with those obtained by filtering the 2 m resolution data over 5 cells in all directions so that the support of the filtered data is also 10 m.Interestingly, in the moderate turbulence case the profiles of the filtered data and unfiltered data are very similar, which is somehow an a posteriori verification that the cutoff lengths associated to the grid-sizes = 2 and = 10 are both situated below the Ozmidov scale in the inertial range of isotropic turbulence.In the weak turbulence, the filtered profiles are still apart and sensibly lower than the corresponding data from the unfiltered data, as the cutoff length of = 10 m is larger than the Ozmidov scale.

Kinetic energy spectra
Given the discretization of physical space, the (discrete) three-dimensional wave numbers in spectral space are where integers n x , varying between 1−N x /2 and N x /2, both inclusive, identify modes number in the direction x, while k x are the corresponding wave numbers: where k x is the wave number spacing.Similar relations hold for y and z directions.The velocity field can be analyzed in spectral space by operating a three-dimensional Fourier transform to the resolved velocity field, where are the coefficients of the Fourier modes that can be calculated by applying a FFT algorithm to Eq. ( 41).Note that in the present approach, periodic boundary conditions are strictly enforced in the two horizontal directions of the computational domain while in the vertical direction, it is the presence of the top and bottom buffer layers that drives the velocity to zero and implicitly insures the periodicity.The energy of the Fourier modes are defined as where the symbol * indicates complex conjugate.The onedimensional spectra of kinetic energy (density of energy per unit wave number) in each of three coordinate directions are obtained by summing up the energy contribution of all Fourier modes that have the same (absolute) mode number in that direction: -5/3 so that for the reality condition, the spectral energy is stored in only half of the modes in each direction.We also define horizontal mode numbers n h , varying between 1 and N h ≡ N x = N y , and vertical mode numbers n v , varying between 1 and N v ≡ N z .The corresponding horizontal and vertical wave numbers (k h , k v ) and wavelengths (λ h , λ v ) are with k h ≡ k x = k y and k v ≡ k z The horizontal energy spectrum is defined as the average of the onedimensional energy spectra in the x and y directions while the vertical spectrum is simply the one-dimensional energy spectrum in the z direction (Lindborg, 2006): Based on dimensional arguments similar to those employed for classical Kolmogorov turbulence, Lindborg (2006) argued that in strongly stratified flows E k (k h ) should scale as k −5/3 h in the inertial range -even though this scaling is associated to the existence of horizontal turbulence cascade rather than to fully developed three-dimensional turbulence.In the inertial range, the theoretical horizontal spectrum then takes the form where C hk is a constant.Figure 11 reports the computed horizontal spectra E k (k h ) for all considered cases along with the corresponding spectra compensated using Eq. ( 51).The range of scales [k hk 1 , k hk 2 ] used to define the inertial range and the scaling factors were determined with a linear regression method and are reported in Table 2.It can be noticed that in the simulations with highest resolutions, = 2 and 4 m, the spectra do show an inertial range approximately constant in the subkilometerscales range, which is of interest to this study, down to the scales on the order of 10 m and below where energy is dissipated by the turbulence model; in particular, for the highestresolution cases, C hk ∼ 0.4 − 0.5, which is on the order of Kolmogorov constant for isotropic turbulence C kol 0.49 (Pope, 2000) starting at scales of around 50 m.The present results are also in line with previous spectral DNS (Lindborg, 2006;Brethouwer et al., 2007;Kimura and Herring, 2012) that considered atmospheric scales larger than the present study.When the forcing is increased, the portion of the spectrum with −5/3 slope tends to increase, which is consistent with the fact that the strength of the turbulence increases relative to the buoyancy force.This tendency can be best appreciated by looking at the compensated spectra in Fig. 11.In the case of weak forcing (and partially for the moderate forcing), when the grid size is decreased to = 10 m, the spectra start to depart from the inertial range at scales much larger than the cutoff length of 2 , which is a marker of excessive dissipation produced by the subgrid-scale model.For stratified flows, this is can be explained by the fact that when turbulence is decreased (or stratification is increased) the thickness of vertical layers also reduces so that finer grid resolution is needed to resolve smaller and smaller structures.In particular, Waite (2011) showed that the buoyancy scale has to be well resolved both in the horizontal and vertical directions.The figure also indicates that the spectra for the highest-resolution cases are slightly shallower than −5/3 with a spectral slope of about −1.45 from scales on the order of 200 m up to 1 km.This tendency was observed for example in DNS by Brethouwer et al. (2007); Kimura and Herring (2012) and Augier et al. (2012) and could be related to the appearance of a spectral bump around the buoyancy wave number (Waite, 2011).To complete the analysis of kinetic energy spectra, Figure 12 shows the evolution of energy in the VSHF modes with k h = 0 (Smith and Waleffe, 2002).In all cases, the energy grows in time at a rate that depends on the resolution and the forcing but in general the tendency is similar to what has been observed in previous DNS (Lindborg, 2006;Smith and Waleffe, 2002).
As for horizontal spectra, it is possible to derive a theoretical expression for the kinetic energy vertical spectra based on dimensional arguments that are valid between buoyancy and the Ozmidov length scales (Billant and Chomaz, 2001;Lindborg, 2006): where C vp is a constant.As noted by Brethouwer et al. (2007), the range of scales where Eq. ( 52) holds is much narrower than the −5/3 inertial range.The vertical spectra E k (k v ) and the compensated spectra E k (k v )/N 2 k −3 v are plotted in Fig. 13 for all the simulated cases.It can be observed that the scaling is pretty well satisfied in the range of wavelengths 50 m < λ v < 200 m.

Potential energy spectra
The spectra of (available) potential energy are obtained from the Fourier transform of the resolved potential temperature deviation from the background state: where the Fourier coefficients are computed using an FFT of θ t .From Eqs. (32), ( 35) and ( 53), the energy of Fourier mode is given by and the horizontal and vertical spectra spectra, E p (k h ) and E p (k v ), are constructed using the same relations as Eqs.( 44)-( 46) and Eqs. ( 49) and ( 50) with the formal substitution of E k by E p .Following similar dimensional arguments used for kinetic energy, the horizontal spectrum of potential energy should scale in the inertial range as (Lindborg, 2006) where C hp is a constant and p is the dissipation rate of potential energy.In the present formulation of the model, the subgrid-scale variance of potential temperature θ 2 is parameterized using Eq. ( 9) rather than being obtained from a dedicated transport equation as it is the case for the subgridscale kinetic energy, Eq. ( 6).Hence, the (local) dissipation p cannot be deduced explicitly from resolved quantities as the dissipation rate of kinetic energy k in Eq. ( 7).The gridaveraged dissipation p is then estimated from the transport equation of potential energy in the Boussinesq approximation, which in statistically steady conditions reads: The horizontal and vertical spectra of potential energy, E p (k h ) and E p (k v ), are reported in Figs. 14 and 15, respectively.The slope of the vertical spectra is slightly shallower than −3 between the buoyancy and Ozmidov scales, which is related to the appearance of the spectral bump as discussed in Sect.3.5.For the horizontal spectra, the figure also shows similar trends observed for the kinetic energy spectra when varying the intensity of forcing except that the inertial range extends over scales larger than the inertial scales of E k (k h ).This difference can be attributed to the fact that potential temperature is not directly forced in these simulations, indeed the depletion of E p (k h ) at very large scales is lower than the corresponding depletion of E k (k h ).The scaling factors C hp of Eq. ( 56) are reported Table 2. Denoting by [k hp 1 , k hp 2 ] the inertial range of potential energy spectra, C hp in this range varies between 0.7 and 0.8 for the high-resolution cases, which is higher than the value C hk ∼0.4-0.5 in Eq. ( 51), observed in the inertial range of kinetic energy [k hk 1 , k hk 2 ].The difference in the location of the inertial range in the kinetic and potential energy spectra can also be attributed to the fact that in the present LES model subgrid-scale energy and temperature are modeled differently, which affects differently the distribution of energy in the resolved scales close to the cutoff.Results from DNS in the literature show that, for a sufficiently high Reynolds number, C hp ∼ C hk over approximately the same inertial range of scales for both kinetic and potential energy spectra (Lindborg, 2006;Brethouwer et al., 2007).In an attempt of comparing the scaling factors in the same way, we calculated the mean value of the compensated spectra of E p (k h ) and E k (k h ) over the same range of scales, i.e., [k hp 1 , k hk 2 ], which covers both the kinetic and potential inertial range.When doing so, the resulting average scaling factors within this range are very close, C hk C hp , and this is for all cases considered.This suggests that, on average, over this range of scales, kinetic and potential energy are transferred down to smaller scales at the same rate although locally the transfer mechanism is different: for kinetic energy it is entirely the result of the turbulent cascade of resolved motions with a constant energy transfer rate (somehow improperly denoted dissipation rate) while for temperature this is due partly to the inertial transfer of resolved motions and partly to the effects of subgrid-scale motions (which act as "viscous" dissipation).To conclude the analysis, the average potential-to-kinetic energy ratio in the range is computed.According to the scaling laws, Eqs. ( 51) and ( 56), E p (k h )/E k (k h ) is exactly equivalent to p / k if the kinetic and potential energy have the same inertial range and scaling factors.In the well-resolved cases, R pk varies between 0.15-0.20 (see Table 2), which is close to p / k ∼ 0.16-0.18,obtained from the values listed in Table 1.The fact that these two ratios are close supports the idea that in the present simulations, the ratio R pk can then be assimilated to the ratio between the potential and kinetic dissipation rates in the range [k hp 1 , k hk 2 ] where C hk C hp .This confirms a posteriori that the choice of the range [k hp 1 , k hk 2 ] is appropriate to represent the average energy transferred in the combined inertial ranges of kinetic and potential energy spectra.

Comparison with observations
To the authors' knowledge there is a lack of observational data of turbulence at the subkilometer scale considered in this study so that a thorough validation is difficult.Despite this, some elements of comparison with observational analysis at -5/3  larger scales are provided in order to verify if the results are in an acceptable range of values.The comparison is mainly based on data sets from the GASP and MOZAIC programs collected from commercial routes.The focus is laid on three aspects of validation: the slope of (horizontal) kinetic energy spectra, the value of turbulence dissipation rate k and (to a limited extent) the variances of turbulent fluctuations.Lindborg (1999) computed the longitudinal and transversal structure functions of the velocity field using MOZAIC data and derived a semi-empirical relation that best fit the data.This analysis showed that for separation distances in the range 10 0 km r 10 2 km, the structure functions scale as r 2/3 , which corresponds to the k −5/3 h scaling of kinetic energy spectra obtained by Nastrom and Gage (1985) using GASP data.Assuming that the spectrum can be extrapolated down to the subkilometer scales, this is in good agreement with the spectra shown in Fig. 11.Some existing observations based on individual flights (Wroblewski et al., 2010), in situ measurements with balloons (Dewan, 1997) and radar measurements (Fukao et al., 1994) also support this scaling in the subkilometer scales.
The average dissipation rate based on all MOZAIC data was k ∼ 6 × 10 −5 m 2 s −3 with some minor variations with respect to the latitude (Lindborg, 1999;Lindborg and Figure 15.Vertical turbulent potential energy spectra (left panels) and compensated spectra (right panels).From top to bottom: strong forcing, moderate forcing, and weak forcing cases.Red lines: = 2 m; green lines: = 4 m; blue lines: = 10 m.

2001
).Similar values were reported by Frehlich and Sharman (2010) from detailed climatology analysis of rawinsonde data and by Fukao et al. (1994) based on routine meteorological observations, and were also reproduced by numerical weather prediction model outputs at scales of 50 km and above (Frehlich and Sharman, 2004).All these data seem to support the results shown in Fig. 4 showing that k = O(10 −5 ) m 2 s −3 .As observed by Lindborg and Cho (2001), these values should not be interpreted as accurate or universal given the intermittent distribution of k in the atmosphere.Indeed, the variance increases when decreasing the averaging length used to calculate from NWP model outputs (Frehlich and Sharman, 2004).Nevertheless, these numbers are useful in that they give an order of magnitude for model comparison.
Finally, the variances of turbulent fluctuations were obtained by Nastrom and Gage (1985) by integrating the energy spectrum for different ranges of wavelengths and averaging for different scenarios (season, latitude, land/see, etc.).As shown, in Table 3, for wavelengths 12.5 km < λ h < 25 km they reported U 2 = 0.388 when averaging all flights, which is in the range of values obtained in the present study.Nastrom and Gage (1985) data are obtained by integrating the portion of energy spectrum in the range of wavelengths 12.5 km < λ h < 25 km, and are reconstructed by summing up the variances in the north-south and east-west directions.

Conclusions and perspectives
High-resolution large-eddy simulations of subkilometerscale turbulence in the upper troposphere-lower stratosphere were carried out with the atmospheric model Méso-NH on computational meshes of up to 8 × 10 9 grid points.Turbulence was sustained by means of a low wave number stochastic forcing method, which allows the flow to reach statistically steady conditions.Turbulence fluctuations and dissipation rates increase with the forcing.In accordance to previous DNS of stably stratified flows, atmospheric stratification leads to a substantially anisotropy of the flow field, which manifests with the presence of elongated horizontal structures.The competition between turbulence and stratification controls the degree of the anisotropy, which increases when forcing is reduced.This was quantified in terms of a local Richardson number, which decreases with the forcing intensity.When forcing is decreased, buoyancy forces tend to overwhelm turbulence, which results in smaller gradients requiring a smaller grid size to be resolved.This impacts the slope of the kinetic and potential energy spectra: for low and moderate forcing, the resolutions of 2 and 4 m are needed to have a correct inertial range with slope close to −5/3 in accordance with theoretical analysis whereas, for strong forcing, 10 m resolution is sufficient.Nevertheless, for strong forcing, the dissipation rate seems too large compared to the DNS results.The vertical energy spectra has a narrow −3 slope between the buoyancy and Ozmidov scales, in accordance with scaling arguments.Considering the scattered and intermittent nature of turbulence in the UTLS, the present results agree reasonably well with available observational analysis.The slope of kinetic energy spectrum, the dissipation rate, and velocity variances are in the range of val-ues obtained from measurements, despite that the available observations are at scales slightly larger than those considered here.This study demonstrated the capability of atmospheric models to reproduce turbulence in the UTLS in the critical subkilometer-scale range.This is also a first step of a more ambitious project that aims at modeling the environmental impact of aviation.In particular, the generated turbulence serves as a background flow-field in follow-up studies addressing the problem of atmospheric dispersion of aircraft emissions and their chemical and microphysical transformations in the free atmosphere, such as the transition of contrails into cirrus clouds.

Figure 1 .
Figure 1.Sketch of the computational domain.

Figure 2 .
Figure 2. Temporal evolution of resolved E kt (upper curves) and E pt (lower curves) for the weak forcing case (red lines), moderate forcing case (green lines), and strong forcing case (blue line).Solid lines: = 2 m; dashed lines: = 4 m; dotted lines: = 10 m.

Figure 3 .
Figure 3. Temporal evolution of the ratio of subgrid-scale to total turbulent kinetic energy R = e k /( E kt + e k ) for the weak forcing case (red lines), moderate forcing (green lines) case, and strong forcing case (blue line).Solid lines: = 2 m; dashed lines: = 4 m; dotted lines: = 10 m.

Figure 6 .
Figure 6.Snapshots of potential temperature fluctuations θ t for the weak forcing case, run S02 (left); moderate forcing case, run M02 (center); and strong forcing case, run S10 (right).The top panels show vertical slides taken at y = L y /2, while the bottom panels show horizontal slides taken at z = L z /2.

Figure 7 .Figure 8 .
Figure 7. Temporal evolution of the buoyancy length L b for the weak forcing cases (red lines), moderate forcing cases (green lines), and strong forcing case (blue line).Solid lines: = 2 m; dashed lines: = 4 m; dotted lines: = 10 m.

Figure 12 .
Figure 12.Temporal evolution of the kinetic energy in the shear modes E k (k h = 0) for the weak forcing cases (red lines), moderate forcing cases (green lines), and strong forcing case (blue line).Solid lines: = 2 m; dashed lines: = 4 m; dotted lines: = 10 m.

Table 1 .
Overview of numerical simulations and principal statistics properties of the flow: run identifier; total number of grid points N grid ; resolution (m); forcing intensity σ f (m s −2 ); turbulent fluctuations σ u , σ v , and σ w (m s −1 ); resolved turbulent potential to kinetic energy ratio E pt / E kt ; subgrid-scale to resolved turbulent kinetic energy ratio e k / E kt ; kinetic energy dissipation rate k (m 2 s −3 ); potential energy dissipation rate p (m 2 s −3 ); buoyancy scale L b (m); Ozmidov scale L O (m), horizontal Froude number Fr.For definitions, seeSects.2.1, 3.2, and 3.4.Data are averaged over the last hour of simulation time during which turbulence is fully sustained by the forcing.

Table 2 .
Definitions of the inertial ranges expressed in wavelengths, λ h = 2 π/k h (m), and scaling factors for kinetic energy spectra[k hk 1 , k hk 2 ]; potential energy spectra [k hp 1 , k hp 2 ]; and combined kinetic and potential energy spectra [k hp 1 , k hk 2 ].The very last column displays the average potential-to-kinetic energy ratio R kp in the range [k hp 1 , k hk 2 ], Eq. (58).N.A. means that no inertial range could be identified.

Table 3 .
Comparison with data from observational analysis.Variances from