# Effects of subgrid-scale turbulence parametrization on the representation of clear-air turbulence using kilometre- to hectometre-scale numerical simulations

Léo Rogel, Didier Ricard, Eric Bazile, Irina Sandu

► To cite this version:

Léo Rogel, Didier Ricard, Eric Bazile, Irina Sandu. Effects of subgrid-scale turbulence parametrization on the representation of clear-air turbulence using kilometre- to hectometre-scale numerical simulations. Quarterly Journal of the Royal Meteorological Society, 2023, 149 (757), pp.3301-3322. 10.1002/qj.4557. meteo-04444266

HAL Id: meteo-04444266

https://meteofrance.hal.science/meteo-04444266v1

Submitted on 23 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Effects of subgrid-scale turbulence parametrization on the representation of clear-air turbulence using kilometre- to hectometre-scale numerical simulations

Léo Rogel<fnref n="1" /> ID | Didier Ricard<fnref n="1" /> ID | Éric Bazile<fnref n="1" /> | Irina Sandu<fnref n="2" /> ID

## Correspondence

Léo Rogel, CNRM (University of Toulouse, Météo-France, CNRS), Toulouse, France

Email: leo.rogel@meteo.fr

A turbulence event arising in a jet exit region above the Belgium–Luxembourg area, determined from airliner in-situ measurements, is reproduced using the meteorological models Application de la Recherche à l'Opérationnel à Méso-Echelle (AROME) and Meso-NH at horizontal resolutions of 1.3 km and 260 m. The behaviour of the subgrid turbulence scheme at 1.3 km and its sensitivity to various parameters are analysed, with results being evaluated using measurements. An increase of the vertical resolution around the tropopause levels with $\Delta z \le 300$ m is shown to greatly enhance the turbulence representation. The use of a non-local formulation of the mixing length in the current parametrization at 1.3 km allows reproduction of a turbulence signal in agreement with the observations. On the contrary, the use of a fully three-dimensional formulation has no impact on the simulation at this resolution (1.3 km). Using the 260 m runs, this turbulence event is linked to hydrodynamical wind shear instabilities characterized by horizontal wavelength of 4.5 km, sub-resolved at the operational resolution. At these small grid-size scales, turbulence evolution and equation budgets reflect an equilibrium between dynamical production and turbulence dissipation, and highlight the importance of horizontal gradients. Subgrid turbulence intensities are assessed to be underestimated by the current parametrization at 1.3 km when compared with this high-resolution reference simulation. Finally, different tests on the turbulence parametrization illustrate a transfer between resolved and subgrid kinetic energy in the model. This transfer stresses the importance of a trade-off between mixing intensity and the representation of wind at resolved scales for the upper troposphere.

## KEYWORDS

clear-air turbulence, high-resolution simulation, shear instabilities, parametrization, upper level jet

# 1 | INTRODUCTION

Clear-air turbulence (CAT) refers to the rapid and chaotic fluctuations of wind in the absence of convective clouds in the free atmosphere (Sharman and Lane, 2016, p. 10).

These turbulent eddy motions can affect aircraft and mostly appear through “bumpiness in flight”. Turbulent events are very active not only in the planetary boundary layer, but also at upper levels near the tropopause, where they reach very high intensities. Sources of CAT at upper

<footnote n="1">¹CNRM (University of Toulouse, Météo-France, CNRS), Toulouse, France</footnote>

<footnote n="2">²European Centre for Medium-Range Weather Forecasts (ECMWF), Reading, UK</footnote>
levels mainly include local wind shear associated with jet stream or upper air fronts, static instability of air layers in the vicinity of convective clouds, and orographic or non-orographic (for instance, convectively induced) gravity wave activity and breaking (Sharman *et al.*, 2012b). These diverse origins can also interact, rendering the dynamics of turbulence events even more complex.

Mid and upper tropospheric CAT forecasts are of crucial interest for commercial airline companies and general aviation, and this interest will likely increase, as the number of CAT occurrences is expected to rise in the context of global warming (Williams, 2017; Prosser *et al.*, 2023). Because CAT represents a serious hazard for aircraft, avoidance strategies based on weather predictions are used to reduce air-passenger injuries, structural damage, and maintenance costs (Sharman and Lane, 2016). Another direct application of CAT forecasts could concern the computation of optic refraction indices from atmospheric turbulence, of particular importance for ground astronomy studies (Masciadri *et al.*, 1999; Masciadri *et al.*, 2017).

Although they have been studied for several decades, CAT phenomena are still particularly difficult to forecast, owing to the range of scales involved and their complex interactions. Early studies investigated the link between CAT and shear instabilities, like Kelvin–Helmholtz billows, that develop in sheared layers of the large-scale atmospheric flow (Reiter, 1969; Dutton and Panofsky, 1970). Conversely, CAT phenomena are known to impact features of the flow at larger scales. For instance, the importance of the role of turbulent fluxes in the evolution of frontal upper level zones and large-scale jets is a result that has been well established for a long time. Early studies used aircraft observations (Kennedy and Shapiro, 1975; Kennedy and Shapiro, 1980) to determine budgets of turbulent quantities and the impact of turbulent fluxes on vertical wind shear and frontogenetical processes. By using research aircraft data, Shapiro (1976) emphasized the primary role of turbulent fluxes in the evolution of potential vorticity in jet systems, especially in the jet exit region, which exhibits enhanced shear due to the curvature of the flow. This analysis was later refined using an isentropic model that parametrized subgrid turbulent fluxes with a K-gradient formulation (Gidel and Shapiro, 1979).

With the increase in computing power over recent decades, current state-of-the-art numerical weather prediction (NWP) systems are now able to explicitly represent synoptic and parts of the mesoscale phenomena, and the related scale interactions. However, because of its microscale nature, upper level turbulent phenomena representation in NWP models still remains challenging at current horizontal and vertical operational resolutions.

One major difficulty in this field is the lack of high-resolution data and measurements describing turbulence at the upper troposphere-lower stratosphere (UTLS) levels. This, to a large extent, is due to the technical difficulties of achieving measurements at such altitudes. On the contrary, in the stable atmospheric boundary layer, numerous studies have investigated the representation of turbulent fluxes (e.g., (Beare et al., 2006; Holtslag et al., 2013; Mahrt, 2014; Couvreux et al., 2020)). Resorting to large-eddy simulation (LES) modelling and experimental campaigns for the characterization of stably stratified turbulent phenomena in the surface layer and above, these studies improved our understanding of the stable regime and enabled substantial enhancements of turbulence parametrization for the boundary layer. However, numerous models use the same parametrization for turbulent fluxes in their upper levels, but these altitudes have not drawn as much attention as the boundary layer has. These altitudes have predominantly strong stable conditions, for which NWP model parametrizations are known to suffer from spurious collapses of turbulence (Canuto et al., 2008, and references therein). Previous results tended to show that the current turbulence parametrization used at Météo-France seems to underestimate turbulent mixing at upper levels (Bazile et al., 2011). It is commonly thought that these deficiencies can be attributed to poorly represented subgrid processes in current parametrizations (Sandu et al., 2013).

For high-altitude CAT events, the predominant source of information that can be used to assess model performance is given by automated in-situ measurements, such as automated eddy dissipation rate (EDR) reports from commercial aircraft, combined with pilot reports, indicating specific occurrences of CAT phenomena. Several case studies of turbulence events, and investigations into their origins, are based on such observations. For instance, Kim and Chun (2010) reported multiple CAT events occurring over the South Korea peninsula and reproduced the meteorological situation with the Weather Research and Forecasting (WRF) Advanced Research WRF model (Skmarock et al., 2008). The authors concluded on the processes conducive to CAT: strong wind shear below the jet core and in the upper flank of the jet, as well as orographic wave breaking above the peninsula. Sharman et al. (2012a) investigated a case of a breaking mountain wave over Greenland and reproduced the incident waves and subsequent breaking signalled by the turbulence parametrization with the Coupled Ocean/Atmosphere Mesoscale Prediction System model at a horizontal resolution of 5 km. Kim et al. (2014) studied near-cloud-turbulence events occurring in an upper level shear layer induced by a precipitation system between Tokyo and Hawaii indicated from pilot reports. Trier et al. (2020) considered
numerous CAT events associated with a large midlatitude cyclone over the United States. Analysed causes include mesoscale gravity-wave breaking with horizontal length scale $\lambda \sim 100$ km, and widespread shear layers, where small-scale perturbations develop into Kelvin–Helmoltz billows with scale $\lambda \sim 5–10$ km. In a companion paper, (Muñoz-Esparza et al., 2020) analysed this case with respect to the current turbulence parametrization in the WRF model running at $\Delta x = 1$ km, comparing it with experiments run at $\Delta x = 370$ m. They found that, depending on the turbulence closure used, wind fluctuations at the UTLS were either too strong or too damped, which had an impact on the representation of upper level mesoscale gravity waves.

Investigations regarding the small-scale features of turbulence are of primary importance for the characterization of CAT events. Thus, studies evaluating the turbulence representation at the UTLS using an NWP model also resort to high-resolution in-situ data from research aircraft. Recent case studies occurring during measurement campaigns include Bramberger et al. (2020), who examined a case above Iceland associated with breaking mountain waves at tropopause levels during the North Atlantic Waveguide and Downstream Impact Experiment. Similarly, Rodriguez Imazio et al. (2022) investigated a CAT event forced by mountain waves above the Drake Passage located in a region of enhanced vertical shear associated with an upper level front, during the Southern Hemisphere Transport, Dynamics, and Chemistry (SOUTHTRAC) campaign. The SOUTHTRAC data have also been used to evaluate new CAT diagnostics systematically from the global Integrated Forecasting System model (Dörnbrack et al., 2022). These observational studies emphasize the particular features of small-scale CAT at upper levels: intermittency with the presence of “bursts” of turbulence, and anisotropy of the wind fluctuations.

Short-range forecasting of CAT has historically been built around the development of diagnostics, computed from large-scale meteorological conditions well resolved by operational models, to localize the presence of turbulence in the area. These operational diagnostics have been designed in parallel with research into the mechanisms leading to turbulence generation. During the last decade, the development of the Graphical Turbulence Guidance (GTG) system at the National Center for Atmospheric Research (Sharman et al., 2006; Sharman and Pearson, 2016) and its operational use has set a milestone regarding turbulence event forecasting. The GTG system uses a weighted aggregation of different user-selected diagnostics on the same metric, and directly updates the weighting method using the available information and derived climatologies from current turbulence reports.

This combination allows one to take into account the numerous sources of turbulence described by the different diagnostics and has been shown to significantly increase forecast accuracy (Sharman and Pearson, 2016). A similar algorithm has been developed at Météo-France, which uses, among other diagnostics, direct outputs from the turbulence parametrization (Crispel and Jaron, 2020), most notably the turbulent EDR index related to the turbulent dissipation from the parametrization: $EDR = e^{1/3}$. The German Weather Service also recently developed a turbulence forecast product, primarily based on the dissipation rate that is computed from the turbulence parametrization scheme (Goecke and Machulskaya, 2021) and supplemented from other quantities that intervene in the turbulence scheme.

The development of new forecasting techniques for CAT thus further motivates improvements of the subgrid turbulence scheme, given its direct impact on these diagnostics. This study is a preliminary step towards this goal, by investigating the behaviour of the operational turbulence scheme in the Météo-France atmospheric models Application de la Recherche à l'Opérationnel Â Meso-Echelle (AROME) and Meso-NH, on the situation of a high-altitude jet above Belgium and Luxembourg. In this article, sensitivity experiments on a specific CAT situation are achieved for both models, with different horizontal resolutions and parametrization settings. The results are qualitatively compared with in-situ EDR measurements and high-resolution simulations to assess the parametrization performance at the operational resolution of 1.3 km.

Section 2 presents the case used throughout the study. A brief analysis of forcing mechanisms for CAT is given, although without going into too much detail, since this is not the main subject of this article. Next, the measurement datasets, models, and parametrizations used are presented. Section 3 first highlights the main results from the sensitivity study on vertical resolution, mixing length formulation, and the choice of one- or three-dimensional turbulence parametrization. Then, comparisons are achieved with Meso-NH runs at 260 m horizontal resolution, in which these high-resolution simulations are used as a reference. Section 4 highlights the main results of the study and presents perspectives for future work.

## 2 | CASE PRESENTATION AND NUMERICAL SET-UP

On January 27, 2018, several in-situ measurements show a moderate-or-greater (MOG; $EDR \ge 0.2 \text{ m}^{2/3} \cdot \text{s}^{-1}$) rated turbulence signal over the Belgium–Luxembourg area between 1200 and 1500 UTC, at around 200 hPa (Figure 1).
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /><figcaption>FIGURE 1 ARPEGE analysis valid at 1200 and 1500 UTC: 250 hPa geopotential (black contours) and wind (shading), and 200 hPa divergence (−6 × 10⁻⁵ s⁻¹, red dashed contour) on January 27, 2018, with available non-null classified eddy dissipation rate in-situ data above 400 hPa (left panel: 1130 to 1230 UTC reports; right panel: 1430 to 1530 UTC reports). Green box: AROME operational, domain 1 (see Table 1 for description); blue box: Meso-NH domain at 1.3 km horizontal resolution, domain 2; red box: Meso-NH domain at 260 m horizontal resolution, domain 3; black box: budget domain. L: low-pressure system; H: high-pressure ridge.</figcaption></figure>

## 2.1 | Characterization of available EDR measurements

EDR in-situ values are derived from equal-duration time series of aircraft vertical acceleration, measured by commercial airliners. Based on Crispel and Jaron (2020), thresholds for this study are as follows: very weak (0.01 ≤ EDR ≤ 0.1 m²/³ ·s⁻¹), weak (0.1 ≤ EDR ≤ 0.2 m²/³ ·s⁻¹), moderate (0.2 ≤ EDR ≤ 0.3 m²/³ ·s⁻¹), and severe (EDR ≥ 0.3 m²/³ ·s⁻¹). The algorithm implemented uses spatial structure functions of the wind components normal to the flight direction, as described in the work of Meischner et al. (2001), who used 100 Hz frequency in-situ data from research aircraft to evaluate Doppler radar detection of thunderstorm turbulence. For the observation data discussed here, measurements come from commercial aircraft. Measurements are taken at a frequency of 1 Hz, and EDR is computed using a rolling average window of 10 s (0.1 Hz). A report is sent every 10 min containing the maximum EDR on the period, as well as the location of the measurement (time, longitude, latitude, and flight level). Moreover, when MOG values are derived from the measurements (EDR ≥ 0.2 m²/³ ·s⁻¹), a report is immediately sent, and then every minute for 6 min.

Although EDR in-situ automated reports allow one to have an objective measurement of turbulence intensity, it may exhibit important biases. Assuming a typical flight speed of 200 m·s⁻¹, these EDR estimates characterize fluctuations at most in the 200 m–2 km scale range. The sampling of these measurements is low, and thus underestimation of EDR should be expected (Wacławczyk et al., 2020). Finally, it is important to recall that commercial aircraft tend to avoid CAT regions when detected.

Available measurements over January 27 include a total of 8,594 automated non-null EDR reports, among which are 7,498 reports over western Europe, between approximately 0000 UTC and 2100 UTC. Figure 2 presents the altitude repartition of the observation dataset used for this study. There are 6,613 reports above 500 hPa altitude, with 179 MOG turbulence intensity reports. Around 76.9% of non-null observations are above 500 hPa, for flights at cruise levels near the UTLS. Beneath 800 hPa there are 2,264 non-null reports with 448 MOG values. On the defined box in Figure 1, there are 11 MOG reports from 1200 UTC to 1615 UTC, above 350 hPa, at an average level of 213 hPa. The altitude histograms (Figure 2) highlight both the boundary layer (800 hPa and below) and the tropopause region (400–150 hPa) as the two main locations of turbulence occurrences. However, contrary
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>FIGURE 2 Altitude repartition of all available non-null eddy dissipation rate reports on January 27, 2018.</figcaption></figure>

to lower levels, aircraft at tropopause flight levels have increased manoeuvring capabilities; for instance, by modifying their cruising altitude. As a consequence, biases linked to avoidance strategies should be more pronounced at these levels.

## 2.2 | Synoptic flow and forcing mechanisms at large scales

At high altitude, a cut-off low-pressure system isolates over the Balearic Islands (L) and then shifts towards Africa during the day. Anticyclonic conditions over France are associated with a powerful ridge (H) over the near Atlantic. A strong upper level jet is present at 250 hPa levels over England which extends over France and Spain by curving under the influence of the upper-level trough (L) over the Mediterranean. Figure 1 shows the 250 hPa wind and geopotential analysis at 1200 UTC and 1500 UTC obtained from the Météo-France global model Action de Recherche Petite Échelle Grande Échelle (ARPEGE) (Courtier et al., 1991). Turbulence manifestations are visible in the region of this upper level jet, revealed by the aforementioned in-situ data from airliners. EDR reports, described in Section 2.1, document two geographic locations experiencing turbulent activity described on the large-scale ARPEGE analysis: one above Belgium, in the prolongation of the Atlantic ridge, down-flow of a strongly anticyclonically curved jet exit region; and a second above the Alps, in the diffuence region of the jet with the presence of strong orography. The study will focus on turbulence over the Belgium area, as this event is signalled by a larger group of well-localized observations.

Several diagnostics were computed from the operational AROME (Seity et al., 2011), corresponding to

the ARO-L90 experiment described in Section 2.3.3. Figure 3a,b shows the two-dimensional adiabatic frontogenesis function F2D (Bluestein, 1992, eq. 2.3.8), and the Lagrangian Rossby number (Zhang et al., 2000) computed as $R_0^L = V^{ageo}/V$, where $V^{ageo}$ is the ageostrophic wind modulus component and $V$ is the wind modulus. Both diagnostics were computed on the 200 hPa isobaric surface. Figure 3c,d shows low values of the gradient Richardson number and negative values of potential vorticity (PV) on the 310 K isentropic surface at 1200 UTC and 1400 UTC. The gradient Richardson number is computed as

$$Ri_g = \frac{N^2}{S^2}, \qquad (1)$$

which represents a local ratio between the Brünt–Vaïsala $N^2$ and vertical shear $S^2$ acting on the fluid parcel. Regions of low $Ri_g$ are indicative of locations where shear instabilities, and consequent turbulence, can arise (Miles, 1961; Miles, 1963; Sharman et al., 2012b). Regions with negative PV are indicative of regions of the flow where inertial instability can develop (Hoskins, 1974).

Based on Figures 1 and 3, upper level convergence and frontogenesis seem to occur upstream of the Belgium region at 200 hPa. This favours the development of sloping stable layers with both enhanced wind shear and reduced static stability $N^2$ on the upper flank of the jet (Ellrod et al., 2003). The strongly curved flow also locally enhances ageostrophic vertical shear, independently of frontogenesis (Knox, 1997). From Figure 3c,d, low $Ri_g$ values below 0.3 illustrate this region of low dynamic stability. Moreover, Figure 3 represents spatial localization of mesoscale waves visualized on negative PV values at 1200 and 1400 UTC. These waves seem to propagate north-eastward, aligned with the ridge direction presented in Figure 1. They develop in the curved anticyclonic region of the jet, where ageostrophy is strong as diagnosed from the Lagrangian Rossby number $R_0^L$. These waves have a corresponding horizontal wavelength of approximately 60–70 km. In such a synoptic configuration, inertia-gravity waves can indeed be spontaneously emitted and propagate away from mesoscale instabilities developing in the unbalanced flow (Plougonven and Zhang, 2014). More importantly, the southern section of the wave structures on the PV propagates through the region of low Richardson symbolized in blue in the figure until 14:00 UTC. This analysis suggests that the CAT event is linked to the propagation of mesoscale inertia-gravity waves into the low $Ri$ region over Belgium. As described in the literature, (Ellrod et al., 2003; Koch et al., 2005), these waves modulate local $Ri$ values, which can thus fall beneath the dynamic instability threshold and allow for the amplification of small-scale shear instabilities.
<figure><img src="image_6.png" /><figcaption>FIGURE 3 Diagnostics computed for ARO-L90 run: (a) two-dimensional frontogenesis function F2D (K·m⁻¹·s⁻¹) and (b) Lagrangian Rossby number $R_0^L$ computed on 200 hPa isobaric surface at 1200 UTC. Negative values of potential vorticity (PV; red), isobars (5 hPa, black contours), and Richardson number $Ri$ (0.3, blue contour) on the 310 K isentropic surface at (c) 1200 UTC and (d) 1400 UTC. Fields have been smoothed using a Gaussian filter of width $5Δx$ for PV and $Ri$ fields, $7Δx$ for the pressure field, and $8Δx$ for $R_0^L$.</figcaption></figure>

This region thus exhibits well-documented large-scale features leading to CAT instabilities, linked to reduced dynamic stability and potential mesoscale gravity wave activity (Uccellini and Koch, 1987; Knox, 1997; Ellrod et al., 2003; Plougonven and Zhang, 2014). Forcing by orographic waves has been inspected using horizontal and vertical sections of vertical wind (not shown), indicating that this mechanism does not seem to play a role for this CAT event.

## 2.3 | Methodology and numerical models

Numerical simulations are performed with the Météo-France operational model AROME (Seity et al.,

2011) and with the mesoscale research model Meso-NH (Lac et al., 2018) at operational grid resolution (1.3 km) for both models and also at 260 m resolution for Meso-NH.

### 2.3.1 | Models

AROME (Seity et al., 2011) is the operational limited-area model of Météo-France, currently running at a $Δx$ = 1.3 km horizontal resolution (Brousseau et al., 2016). AROME solves a fully compressible equations set (Bubnová et al., 1995) over a spectral grid. The model is integrated using a semi-Lagrangian, semi-implicit method for all prognostic variables, including turbulent
<figure><img src="image_7.png" /></figure>

kinetic energy (TKE) with an operational time step of $\Delta t = 50$ s. The vertical discretization follows a hybrid pressure terrain-following coordinate (Laprise, 1992). Lat-eral boundaries are enforced by the ARPEGE large-scale model. The numerical integration also uses a spectral trun-cation of orography at $3\Delta x$, resulting in a loss of variance that can be distinguished on the tail of kinetic energy spectra (e.g., Ricard et al., 2013).

Meso-NH (Lac et al., 2018) is a non-hydrostatic mesoscale research model jointly developed by the Labo-ratoire d'Aérologie and the Centre National de Recherches Météorologiques. In this study, Meso-NH is configured to solve the pseudo-incompressible system of Durran (1989). The model is integrated using an explicit fourth-order Runge-Kutta scheme in time and centred fourth-order advection scheme for the momentum equations. All scalar quantities, including the TKE, are advected using a piecewise parabolic method scheme (Colella and Wood-ward, 1984), and the vertical discretization follows the orography (Gal-Chen and Somerville, 1975). For all Meso-NH runs, the time step is set to $\Delta t = 2$ s.

Both AROME and Meso-NH use the same physics, including microphysics parametrization (ICE3 scheme; Pinty and Jabouille (1998)), shallow convection parametrization (PMMC09 scheme; Pergaud et al. (2009)), turbulence parametrization (CBR scheme; Cuxart et al. (2000)), radiation scheme (Fouquart and Bon-nel, 1980; Mlawer et al., 1997), allowing the attribution of discrepancies between their results to the different dynamics and numerical methods when they are used in the mesoscale configuration. The AROME operational suite notably uses an enhanced numerical diffusion on the vertical (evolving as $1/p$ on the vertical, with $p$ the pressure) and a much higher time step than the research model. This leads to more diffusive numer-ics in the operational model and a less fine effective resolution when compared with the research model (Ricard et al., 2013).

Finally, all models use the SURFEX module for sur-face coupling. The orography is obtained from the 250 m resolution SRTM and the ECOCLIMAP II v.2.3 databases. Data are adapted according to the resolution of the model used. Surface fluxes are computed using the different schemes ISBA (Noilhan and Planton, 1989), TEB (Mas-son, 2000), ECUME (Belamari and Pirani, 2007), and FLAKE for Meso-NH runs, according to the type of cover (respectively continental surface, towns, sea, and lakes). These parametrizations are expected to have an impact limited in the lowest levels of the models. Moreover, the subgrid-scale (SGS) orography does not contribute to the computation of turbulent fluxes in the atmospheric parametrization.

## 2.3.2 | Turbulence scheme

The models implement an eddy-diffusivity-based turbu-lence parametrization (Cuxart et al., 2000), consisting of a 1.5-order closure scheme based on a prognostic TKE evolu-tion equation, with a diagnostic computation of the mixing length $L$. In the following, subgrid TKE, or SGS TKE, refers to the parametrized TKE computed at the numerical resolution of the model.

In a one-dimensional, single-column framework (hereafter T1D), the subgrid TKE $e$ satisfies the following equation:

$$ \frac{De}{Dt} = \underbrace{\frac{\partial}{\partial z} \left( C_{TT} L \sqrt{e} \frac{\partial e}{\partial z} \right)}_{\text{DT}} - \underbrace{\overline{u'w'} \frac{\partial \overline{u}}{\partial z} - \overline{v'w'} \frac{\partial \overline{v}}{\partial z}}_{\text{DP}} + \underbrace{\beta \overline{w'\theta_v'}}_{\text{TP}} - \underbrace{C_e \frac{e^{3/2}}{L}}_{\text{DISS}}, \quad (2) $$

with $C_{TT} = 0.4$ and $C_e = 0.85$, and with DT being the turbulent transport (acting as a diffusion term), DP the dynamical production of TKE from wind shear, TP the thermal production (positive in convective cases and neg-ative in stable cases), and DISS the turbulence dissipation rate. $u, v$, and $w$ represent the wind components, $\theta_v$ the virtual potential temperature, and $\beta = g/\theta_{vref}$ with $g$ the gravitational acceleration and $\theta_{vref}$ the reference potential temperature profile. Primed quantities denote the fluctua-tions from the ensemble averaged quantities denoted by an overbar. A fully three-dimensional scheme is also imple-mented in Meso-NH (T3D), which integrates the horizon-tal gradients of wind into the computation of the TKE dynamical production.

The operationally used mixing length is the one from Bougault and Lacarrere (1989) (BL89), computed from the upward and downward paths of an air parcel undergo-ing buoyancy and subgrid kinetic energy effects, for both stable and unstable cases,

$$ \int_{z}^{z+l_{up}} \beta(z)(\theta_{v}(z') - \theta_{v}(z)) dz' = e(z), \quad (3) $$

$$ \int_{z-l_{down}}^{z} \beta(z)(\theta_{v}(z) - \theta_{v}(z')) dz' = e(z) \quad (4) $$

and evaluates the master length scale as an average of both $l_{up}$ and $l_{down}$:

$$ L_{BL89} = \left( \frac{l_{up}^{-2/3} + l_{down}^{-2/3}}{2} \right)^{-3/2} . \quad (5) $$
Experiments are also done with the Rodier et al. (2017) (RM17) mixing length, which adds a vertical shear-dependent term S in the integrand of the $l_{up}$ and $l_{down}$ intermediary length scales (Equation (4)):

$$ \int_{z}^{z+l_{\text{down}}} [\beta(z)(\theta_{v}(z') - \theta_{v}(z)) + C_{0}\sqrt{eS}] dz' = e(z), \quad (6) $$

where

$$ S^2 = \left( \frac{\partial \bar{u}}{\partial z} \right)^2 + \left( \frac{\partial \bar{v}}{\partial z} \right)^2. $$

This modification was initially motivated by an overestimation of the TKE in the weakly stable regime, with the integrals becoming null for the neutral regime and thus inducing ill-defined $l_{up}$ and $l_{down}$. In the case of a neutral, constant shear flow, the RM17 formulation is bounded by the Hunt (1988) formulation:

$$ l_H = \frac{1}{C_0} \frac{\sqrt{e}}{S}. $$

Finally, another mixing length specification used is from Deardorff (1980), based on a minimization between the cell size and the local static stability of the flow $N^2$ ($N$ being the Brunt–Väisälä frequency):

$$ L_{\text{DEAR}} = \min \left( (\Delta x \Delta y \Delta z)^{1/3}, 0.76 \sqrt{\frac{e}{N^2}} \right), \quad (7) $$

with

$$ N^2 = \frac{g}{\theta_{\text{vref}}} \frac{\partial \bar{\theta}}{\partial z}, $$

which has been notably used in LESs of stably stratified flows in the boundary layer (Couvreux et al., 2020).

Once the TKE and the mixing length have been updated, turbulent fluxes for the $\phi$ quantity are then computed following the general formulation (for stable cases):

$$ \overline{u_i' \phi'} = -K_{\phi} L \sqrt{e} \frac{\partial \bar{\phi}}{\partial x_i} f_i, \quad (8) $$

where $K_{\phi}$ is a positive constant parameter and $f_i$ a function of stability, dependent on the Redelsperger numbers (Redelsperger and Sommeria, 1981) and only active on the vertical.

### 2.3.3 | Summary of numerical experiments

A series of simulations is used to assess the sensitivity to various parameters at 1.3 km: the mixing length formulation, the use of the three-dimensional scheme for the TKE equation, an increased vertical resolution of the model at tropopause levels, and the reduction of the turbulence dissipation rate. Additionally, high-resolution numerical experiments at $\Delta x = 260$ m are also conducted and compared with the control run (MNH-CTRL) at $\Delta x = 1.3$ km.

The simulations carried out are summarized in Table 1. The ARO-L90 and MNH-L90 runs use the vertical levels of the operational AROME model, with a stretched grid on the vertical and a mean grid size of $\Delta z \approx 400$ m in the $z = 10–12$ km range. The ARO-L105 uses the ARPEGE model 105 levels, which feature an increased vertical resolution at the tropopause layer of $\Delta z \approx 290$ m. The other Meso-NH runs use the same vertical grid

<table><caption>TABLE 1 List of the numerical experiments and their corresponding domain, grid spacing $\Delta x$, average vertical grid spacing $\Delta z$ for $z = 10–12$ km, mixing length formulation (see Section 2.3.2 for definitions), scheme formulation, and $C_{\epsilon}$ value.</caption><thead><tr><th>Label</th><th>Model</th><th>Domain</th><th>Δx (km)</th><th>Δz (m)</th><th>Vertical level numbers</th><th>Mixing length</th><th>Scheme</th><th>C<sub>ε</sub></th></tr></thead><tbody><tr><td>ARO-L90</td><td>AROME</td><td>1</td><td>1.3</td><td>~400</td><td>90</td><td>BL89</td><td>T1D</td><td>0.84</td></tr><tr><td>MNH-L90</td><td>Meso-NH</td><td>2</td><td>1.3</td><td>~400</td><td>90</td><td>BL89</td><td>T1D</td><td>0.84</td></tr><tr><td>ARO-L105</td><td>AROME</td><td>1</td><td>1.3</td><td>~290</td><td>105</td><td>BL89</td><td>T1D</td><td>0.84</td></tr><tr><td>MNH-CTRL</td><td>Meso-NH</td><td>2</td><td>1.3</td><td>~290</td><td>105</td><td>BL89</td><td>T1D</td><td>0.84</td></tr><tr><td>MNH-BL3D</td><td>Meso-NH</td><td>2</td><td>1.3</td><td>~290</td><td>105</td><td>BL89</td><td>T3D</td><td>0.84</td></tr><tr><td>MNH-DE3D</td><td>Meso-NH</td><td>2</td><td>1.3</td><td>~290</td><td>105</td><td>Deardorff</td><td>T3D</td><td>0.84</td></tr><tr><td>MNH-DE1D</td><td>Meso-NH</td><td>2</td><td>1.3</td><td>~290</td><td>105</td><td>Deardorff</td><td>T1D</td><td>0.84</td></tr><tr><td>MNH-RM17</td><td>Meso-NH</td><td>2</td><td>1.3</td><td>~290</td><td>105</td><td>RM17</td><td>T1D</td><td>0.34</td></tr><tr><td>MNH-DISS</td><td>Meso-NH</td><td>2</td><td>1.3</td><td>~290</td><td>105</td><td>BL89</td><td>T1D</td><td>0.42</td></tr><tr><td>MNH-HR1D</td><td>Meso-NH</td><td>3</td><td>0.26</td><td>260</td><td>96</td><td>Deardorff</td><td>T1D</td><td>0.84</td></tr><tr><td>MNH-HR3D</td><td>Meso-NH</td><td>3</td><td>0.26</td><td>260</td><td>96</td><td>Deardorff</td><td>T3D</td><td>0.84</td></tr></tbody></table>
spacings as ARO-L105. Corresponding to the ARO-L105 run, a specific Meso-NH run is defined as a control run (MNH-CTRL). This run uses the same physical parametrization as the current operational AROME model (one-dimensional scheme, BL89 mixing length). For the $\Delta x = 260$ m grid-size experiments MNH-HR1D and MNH-HR3D, the vertical grid size of the model has been constantly set to 260 m between 7.5 and 14 km altitudes in order to enhance the isotropy of the grid at levels of interest.

# 3 | RESULTS

The results are divided according to the resolution at which models were run. In Section 3.1, the sensitivity of the models to various parameters at the operational resolution ($\Delta x = 1.3$ km) is shown. Section 3.2 presents the results at $\Delta x = 260$ m and also compares them with the runs at $\Delta x = 1.3$ km.

## 3.1 | Operational resolution numerical experiments

Figure 4 shows TKE horizontal sections at 225 hPa for the numerical experiments ARO-L105 and MNH-CTRL, with 250 hPa horizontal wind barbs. TKE and wind fields were linearly interpolated on isobaric surfaces for comparisons on the same coordinate levels. Wind velocities in the divergence zone of interest reach $55 \text{ m} \cdot \text{s}^{-1}$. At 1200 UTC, a region of moderate TKE intensity ($\text{TKE} \ge 1 \text{ m}^2 \cdot \text{s}^{-2}$) develops above eastern Belgium, on the northern part of the box in both the ARO-L105 and MNH-CTRL runs, which is in good agreement with the spatial localization of the in-situ reports. Also worth noticing is the higher effective resolution for the Meso-NH model than for AROME (Ricard et al., 2013; Skamarock et al., 2014), the former representing smaller scale patterns in the resolved TKE field.

A more thorough comparison between both simulations and the observations can be achieved using the in-situ measurements. However, because numerical models express $\epsilon$ as a function of TKE, comparison between model EDR values and observations is not direct. The turbulence scheme in both AROME and Meso-NH outputs instantaneous values of TKE at the grid size of 1.3 km. This scale approximately corresponds to the larger scales upon which aircraft are sensitive to turbulent fluctuations. Indeed, aircraft are mostly sensitive to fluctuations on the 10 m–1 km range (Sharman et al., 2014). It is important to be mindful that scales represented by model EDR values and EDR observations can differ in regard

<figure><img src="image_8.png" /><figcaption>**FIGURE 4** The 225 hPa turbulent kinetic energy (TKE) field (shading) and 250 hPa wind barbs ($\text{m} \cdot \text{s}^{-1}$) at 1200 UTC for (a) ARO-L105 and (b) MNH-CTRL runs. The blue line corresponds to the oblique sections in Figure 7.</figcaption></figure>

to the spatial scales they are inferred from. Here, in-situ measurements are only used qualitatively to examine the sensitivity of vertical resolution modifications for the models with respect to observations. As additional information, model-derived values for EDR = $\epsilon^{1/3}$ have been used in other studies to compare kilometre-scale simulation outputs from the subgrid parametrization and measurements inside the boundary layer; for instance, Muñoz-Esparza et al. (2018) compares $\epsilon^{1/3}$ with sonic anemometer data.

In this study, model values for EDR are retrieved using the gridded data of instantaneous values of the dissipation rate $\epsilon = C_{\epsilon}e^{3/2}/L$. These values are obtained after the resolution of the prognostic TKE field from Equation (2).
<figure><img src="image_9.png" /><figcaption>**FIGURE 5** Comparison of eddy dissipation rate $EDR = \epsilon^{1/3}$ from in-situ measurements and model values inside the budget box, between 350 hPa and 150 hPa, from 1000 to 1400 UTC. EDR values below 0.01 $m^{2/3} \cdot s^{-2}$ are filtered out.</figcaption></figure>

Each grid value of the dissipation is then projected on the $EDR = \epsilon^{1/3}$ metric. Finally, a statistical comparison is made with the observations. Model data are classified according to the same thresholds presented on Section 2.1. Since EDR in-situ reports truncate the occurrences when no fluctuations are detected, values obtained by the models were filtered in respect to local TKE values. Here, $TKE \ge 1 \times 10^{-4} m^2 \cdot s^{-2}$ was used. The same filtering is applied for the three experiments ARO-L90, ARO-L105, and MNH-CTRL.

The comparison between model EDR values and measured values is shown in Figure 5, using all available reports between 350 and 150 hPa from 1000 to 1400 UTC over the black box. More than 80% of the data correspond to weak and very weak turbulence measurements, and 18% of the data indicate moderate turbulence. The proportion between weak/very weak and moderate/severe turbulence is fairly well represented by the simulations.

The ARO-L90 run, using the 90 operational levels, does not reproduce the data repartition of the available measurements well. However, switching to the 105 levels of the ARPEGE model (ARO-L105) drastically enhances the model behaviour in that regard. The MNH-CTRL and ARO-L105 runs behave similarly, with around 53% of very weak turbulence, 28.5% of weak, 13% of moderate, and 3% of severe values. ARO-L105 also overestimates the proportion of very weak turbulence and underestimates the proportion of weak and moderate CAT intensity. With the same trend, the MNH-CTRL run better approaches the observations. Both ARO-L105 and MNH-CTRL runs also forecast severe values of CAT that are not present in the observations, and it is not clear if this is due to an overestimation from the models or a lack of representativeness of this CAT event from the in-situ measurements induced

by aircraft avoidance strategies and the small number of observations.

Figure 6 also shows a qualitative comparison of the time evolution of observed (EDR mesurements) and simulated (mean TKE and EDR instantaneous values) turbulence over the selected domain (black box). Observations indicate MOG turbulence between 1030 and 1630 UTC with a peak around 1230 UTC and another around 1500 UTC. The two runs both reproduce an increase of turbulence present in the observations, starting at 1030 UTC. For MNH-CTRL, TKE and EDR both reach peak values at around 1330 UTC. For the ARO-L105 run, TKE is maximum around 1500 UTC and EDR peaks at 1600 UTC.

Spatial vertical localization of TKE with respect to the jet can be determined from the oblique sections in Figure 7 for the MNH-CTRL run. We also represent the local $Ri_g$ number, indicative of the local destabilization of the flow with respect to shear instabilities (Miles, 1961). Most of the TKE is produced in regions of low gradient Richardson number ($\le 0.25$) above the jet at altitudes between 10 and 12 km. Low $Ri_g$ values indicate the unstable feature of the shear layers from the background flow above the jet, which might be conducive to shear instability and turbulence.

Figure 8 shows the vertical profiles of the mean budget terms of the subgrid TKE equation over 1 hr for the MNH-CTRL and ARO-L105 runs, as well as the vertical profile of mean hourly variation of TKE. Since both runs use the one-dimensional formulation of the turbulence scheme, only vertical gradients are used for the integration of the TKE equation. For the hours considered, the TKE time evolution is positive, indicating a TKE production. For both models, the only positive term is the shear-induced dynamical production. The leading sink term is the turbulence dissipation $\epsilon$. The thermal production terms are negative, characteristic of the strongly stable features of the background flow. Finally, the advection and subgrid turbulent transport terms are negligible with regard to the other terms. Thus, in the model, the subgrid TKE is only produced by the jet vertical wind shear. As will be discussed in Section 3.1.2 (and Figure 6), horizontal gradients play no significant role in the TKE production at this resolution.

### 3.1.1 | Sensitivity to vertical resolution

Here, we discuss the sensitivity of the TKE production relative to the vertical grid sizes in the AROME model, and the impact of switching from the default 90 operational levels to the 105 levels of the ARPEGE model. Although the ARO-L105 run manages to produce TKE from 1100 UTC, the beginning of the production in ARO-L90 is delayed up to 1400 UTC (see Figure 6) with a smaller peak around
<figure><img src="image_10.png" /><figcaption>**FIGURE 6** Time evolution in the black box (Figure 1), between 8 km and 14 km of (a) filtered eddy dissipation rate (EDR) reports (the solid curve is a Bezier fit of the points), (b) model EDR (see text for details), and (c) averaged subgrid turbulent kinetic energy (TKE), weighted accordingly to the layer thicknesses of corresponding levels in the model, for all numerical experiments summarized in Table 1.</figcaption></figure>

1600 UTC. As a consequence, this run performs very poorly when compared with observations between 1000 and 1400 UTC (Figure 5), only indicating very weak turbulence.

Figure 9 compares stability values $Ri_g$, $N^2$, and $S^2$ derived from mean wind and temperature profiles for the AROME runs, and the impact of the vertical resolution increase on stability of the resolved flow. With ARO-L105, slightly increasing the number of levels below the tropopause layer changes the vertical shear profile, increasing $S^2$ values between $z = 9,500$ m and 11,500 m. The $\Delta z$ change also impacts the $N^2$ values indirectly through mixing, but the net effect on the $Ri_g$ profile is the mean reduction of stability in the region of interest. This suggests that adding several levels in the UTLS improves the TKE production following the increase of resolution in the vertical gradients. It also highlights that these gradients seem sub-resolved in the operational AROME 90-levels vertical grid. A similar trend is observed with the MNH-L90 run, to be compared with MNH-CTRL. It also shows a delay in TKE production, but to a lesser extent than for AROME runs.

### 3.1.2 | Sensitivity to turbulence scheme parameters

As indicated earlier, other runs used different mixing-length formulations, the dimension of the

turbulence scheme (T1D versus three-dimensional turbulence [T3D]), and the parameter value for the dissipation. The time evolution of the mean TKE and of model EDR for all the experiments at $\Delta x = 1.3$ km is shown in Figure 6. In this figure, model EDR has not been filtered with regard to TKE values, contrary to the results presented in Figure 5.

The influence of vertical gradients that was highlighted in Section 3.1.1 also raises questions regarding that of the horizontal gradients in the TKE computation. The evolution for the MNH-BL3D run (Figure 6), which uses the T3D scheme, clearly shows that, for $\Delta x = 1.3$ km in this case, the full three-dimensional computation for the turbulence scheme does not modify TKE and EDR intensities at the UTLS levels. This is in contrast to several other studies on other parts of the atmosphere; for instance, Goger *et al.* (2018), which highlights the impact of horizontal gradients on the production of TKE in an Alpine valley boundary layer at $\Delta x = 1.1$ km with steep orography, or in the free troposphere, where Verrelle *et al.* (2015) also showed that the difference between T1D and T3D becomes perceptible at 2 km resolution for simulations of deep convection and pointed out the necessity to deal with horizontal turbulent fluxes at kilometre resolutions.

The use of the Deardorff mixing length has a huge impact (Figure 6), with almost no TKE production at all during the entire simulation, and an additional sensitivity test (not shown here) shows that it is the threshold effect of the formulation that inhibits turbulence production. On the EDR time evolution, very low values of dissipation
<figure><img src="image_11.png" /></figure>

**FIGURE 7** Vertical cross-section "A" along the blue line shown in Figure 4 for the MNH-CTRL run of (a) turbulent kinetic energy (TKE) and (b) gradient Richardson number $Ri_g$ number. The temperature field is represented as grey contour lines (3 K isotherms), and the jet core is represented by the blue line ($45 \text{ m} \cdot \text{s}^{-1}$ isotach).

are also retrieved. Also, the effect of the RM17 formulation is to reduce TKE, which is not surprising since it was designed to reduce turbulent mixing in weakly stable cases in the boundary layer. Finally, the experiment MNH-DISS with a reduced turbulent dissipation coefficient in the scheme enhances TKE and mixing. Also for the MNH-DISS run, EDR is enhanced compared with MNH-CTRL before the peak at 1330 UTC, but after the peak it stays below the control run values.

Varying intensities of subgrid mixing also have a direct impact on the representation of gradients and winds in the model. Figure 10 shows several vertical profiles on the same box from the MNH-CTRL, MNH-DE1D, and MNH-DISS experiments. As noticed before, the Deardorff mixing length produces no subgrid TKE and serves as a reference with no turbulent mixing, to be compared with MNH-CTRL and MNH-DISS runs. The mean TKE in the box is higher for MNH-DISS, which also produces values in the lower flanks of the jet, between 7.5 and

<figure><img src="image_12.png" /></figure>

**FIGURE 8** Vertical mean profiles of (a) budget terms of the turbulent kinetic energy (TKE) equation between 1100 and 1200 UTC over the black box shown in Figure 3 and (b) TKE difference between 1100 and 1200 UTC over the same box. ADV: advection; DP: dynamical production of TKE from wind shear; DISS: turbulence dissipation rate; DT: turbulent transport; TP: thermal production.

9.5 km. Mean mixing length reaches values up to ~ 79 m for MNH-CTRL and ~ 108 m for MNH-DISS. Another feature of the enhanced subgrid mixing is the overall reduction of gradient intensity, especially between 10 and 13 km, where TKE is the strongest. At these levels, $Ri_g$ is higher for MNH-DISS and lower for MNH-DE1D, with the highest impact coming from the reduction of $S^2$ in the former. On the wind profiles, the smoothing effect of increased turbulent mixing can also be noticed for both zonal and longitudinal components. This reflects a balance between the resolved part and subgrid part of the kinetic energy in the model.

## 3.2 | High-resolution numerical experiments

The following examines results from the $\Delta x = 260 \text{ m}$ runs with Meso-NH. They are also compared with the results from the MNH-CTRL run at operational resolution.

### 3.2.1 | Main features of CAT at high resolution

At high resolution, the model resolves small-scale wave-like structures visible in TKE, and vertical wind
<figure><img src="image_13.png" /><figcaption>**FIGURE 10** Comparison of (a) gradient Richardson number $Ri_g$, (b) Brünt–Vaïsala frequency $N^2$ (c) vertical shear $S^2$, (d) turbulent kinetic energy (TKE), (e) mixing length, and (f) wind components (dashed lines for $U$, solid lines for $V$) for MNH-CTRL, MNH-DE1D, and MNH-DISS runs at 1200 UTC.</figcaption></figure>

<figure><img src="image_14.png" /></figure>

fields, as can be seen in Figure 11, which are orthogonal to the direction of the jet streak. These structures have a corresponding horizontal wavelength of 4.5 km, which also corresponds to wavelengths of wind shear

instabilities, as will be discussed further later. Based on 10 min outputs of vertical winds, structures appear in the zone of interest (black box) around 1100 UTC and evolve eastward with the synoptic flow to exit this area at
<figure><img src="image_15.png" /><figcaption>FIGURE 11 (a) The 250 hPa turbulent kinetic energy (TKE; red shading) and 225 hPa wind barbs for MNH-HR3D; (b) same, but with 250 hPa vertical wind (WT; blue-red shading).</figcaption></figure>

around 1500 UTC. With the use of oblique cross-sections (Figure 12), we can spot these perturbations on the same levels of reduced $Ri_g$ as for the experiments at coarser resolutions, and perturbations of $\theta$ extending up to around 18 km in the stratosphere. Locally, the vertical wind velocities can reach up to $7-8 \text{ m}\cdot\text{s}^{-1}$ at 200 hPa. These strong values are a particularly striking feature for these CAT events, occurring outside of convection and in regions where classical large-scale assumptions are not valid. In contrast, values for large-scale vertical wind velocity estimations in geostrophic equilibrium are of the order of $w \sim 1 \text{ cm}\cdot\text{s}^{-1}$ (Bluestein, 1992, section 5.2.2). The high values of vertical wind obtained here are analogous in magnitude to the strong values obtained for CAT events from in-situ measurements (Kennedy and Shapiro, 1980; Bramberger et al., 2020), which reported values in the range of $w \sim 5 \text{ m}\cdot\text{s}^{-1}$. Moreover, recent numerical studies of CAT events (Trier et al., 2020), as well as more idealized frameworks (Mahalov et al., 2007), have reported a similar

<figure><img src="image_16.png" /><figcaption>FIGURE 12 As in Figure 7 for MNH-HR3D at 1330 UTC, but for cross-section “B” in Figure 11.</figcaption></figure>

range of values for vertical velocities generated by wind shear instabilities in the jet flanks.

The sensitivity of the scheme to horizontal wind gradients at this resolution is highlighted in Figure 13, where mean budget terms during half an hour at different instants are shown, at time locations indicated in Figure 14 with grey vertical bands. We can notice that the T3D formulation of the scheme (solid lines) nearly doubles TKE intensities (Figure 14) and the dynamical production term of TKE compared with T1D (dashed lines). In the early stage of TKE production (e.g., see between 1300 and 1330 UTC), dynamical production is slightly larger than dissipation, allowing the increase of subgrid turbulence. At subsequent times, the subgrid TKE evolution is mostly governed by an equilibrium between dynamical production and turbulent dissipation, as the other terms remain smaller. In particular, the thermal production is negligible in the core production layer $z = 11-12$ km, where

$$ \left| \frac{\text{TP}}{\text{DISS}} \right|_{\text{max}} \sim 0.05, \qquad (9) $$
<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /><figcaption>FIGURE 13 Mean budget terms of the turbulent kinetic energy (TKE) equation for MNH-HR1D and MNH-HR3D run on the budget box between (a) 1300 and 1330 UTC, (b) 1330 and 14:00 UTC, and (c) 1430 and 1500 UTC. ADV: advection; DP: dynamical production of TKE from wind shear; DISS: turbulence dissipation rate; TP: thermal production; TR: turbulent transport.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>FIGURE 14 Comparison of the coarse-grained and subgrid turbulent kinetic energy (TKE) at 1.3 km horizontal scale from the MNH-HR1D and MNH-HR3D runs with MNH-CTRL run, between 8 and 14 km, on the same budget domain. SGS: subgrid-scale.</figcaption></figure>

contrary to the $\Delta x = 1.3$ km budgets in Figure 8 for MNH-CTRL, for which

$$ \left| \frac{\text{TP}}{\text{DISS}} \right|_{\text{max}} \sim 0.25. \qquad (10) $$

This suggests that the subgrid part of the turbulence due to thermal production decreases with higher resolution. This is in line with results obtained for deep convection showing that the subgrid thermal production decreases with increasing resolution, whereas the part of the subgrid dynamic production increases (Strauss et al., 2022).

### 3.2.2 | Comparison of small-scale variability across simulations at different resolutions

To compare experiments across different horizontal resolutions, it is preferable to take one scale as a reference, since subgrid TKE does not represent the same amount of kinetic energy depending on the scale at which it is computed. At the 1.3 km resolution, the subgrid TKE is assumed to represent kinetic energy that is below this scale. With this perspective in mind, we use a coarse graining method (Honnert et al., 2011; Dorrestijn et al., 2013), which consists of aggregating on coarser horizontal boxes the high-resolution datasets. This defines an averaging operator that is used to compute fluctuations and variances from the mean values in each of the boxes. The subgrid values of the TKE in the high-resolution dataset must also be taken into account, since it is also subgrid at the coarser resolution. Denoting the coarse-grained (CG) quantity $\phi$ as $\overline{\phi}^{\Delta x}$, the CG TKE $e$ at an arbitrary horizontal resolution $\Delta x$ is computed as

$$ \underbrace{e(\Delta x)}_{\text{CG TKE}} = \underbrace{\overline{e_{\text{HR}}^{\Delta x}}}_{\text{SGS TKE}} + \frac{1}{2} \sum_{i=1,3} \overline{(u_i - \overline{u_i}^{\Delta x})(u_i - \overline{u_i}^{\Delta x})}^{\Delta x}, \quad (11) $$

with $e_{HR}$ being the parametrized TKE from the high-resolution grid on MNH-HR1D and MNH-HR3D runs. Figure 14 presents the temporal evolution of the CG TKE at $\Delta x = 1.3$ km as well as the SGS part of this diagnostic on the box, also between 8 and 14 km for the
MNH-HR1D and MNH-HR3D experiments. For the sake of comparison, the TKE from the MNH-CTRL run is also presented. This diagnostic indicates that at 260 m resolution the resolved wind field contains more horizontal variability than at 1.3 km. As noted before, the T1D SGS TKE is nearly half that of the T3D run at $\Delta x = 260$ m, and this has repercussions on the diagnosed CG TKE from MNH-HR1D and MNH-HR3D, the latter being higher in intensity after the peak production of SGS TKE at 1330 UTC. The CG TKE of the two runs experiences a rapid growth, compatible with an exponential curve of the form $e(\Delta x = 1.3 \text{ km}) \propto \exp(-\omega t)$, with $\omega$ being the growth rate of the perturbations. Values of this parameter are discussed later. When compared with the diagnosed CG TKE, MNH-CTRL underestimates turbulent mixing at the UTLS.

The discrepancies between resolutions can also be seen in the horizontal power spectral density of wind obtained from a discrete cosine transform (DCT) (Denis et al., 2002; Ricard et al., 2013). Denoting $\hat{\phi}^{\text{DCT}}$ the DCT transform for a quantity $\phi$, we note $E_h = \hat{U}^{\text{DCT}} + \hat{V}^{\text{DCT}}$ the spectral density for horizontal winds and $E_z = \hat{W}^{\text{DCT}}$ for vertical winds. Figure 15 presents $E_h$ and vertical winds $E_z$ for various experiments at 1.3 km and at 260 m grid resolutions, and averaged between the 175, 200, and 225 hPa levels. A peak of kinetic energy can be distinguished on the computed spectra for runs MNH-HR1D and MNH-HR3D at the horizontal scale of 4.5 km. In scales smaller than this energy peak, the development of small-scale structures is noticeable that appear to undergo a turbulent energy cascade following a $\kappa^{-5/3}$ law, with $\kappa$ being the horizontal wave number. According to Kolmogorov (1941), in the presence of isotropic, homogeneous turbulence, the three-dimensional spectrum of small-scale motions should follow the Kolmogorov law $E(\kappa) \sim K\epsilon^{1/3}\kappa^{-5/3}$ with $\epsilon$ constant across the scales in the inertial range and $K = 1.5$ a dimensional constant. However, for stable conditions around the tropopause levels, one must rather assume that turbulent fluctuations are anisotropic, which brings questions concerning the kinetic energy cascade process. This is discussed further in Section 3.3. Also worth noting is the fact that the MNH-HR3D run exhibits less energy content in the small scales, supposedly due to the increased mixing from the turbulence scheme when accounting for horizontal wind gradients.

On the contrary, most of the runs at the $\Delta x = 1.3$ km do not exhibit this high variability on the smaller resolved scales. Only the MNH-CTRL and the MNH-DE1D runs seem to partially resolve the energy peak at 4.5 km. The MNH-CTRL splits from the high-resolution spectrum at $\lambda \sim 13$ km for the vertical wind and at 20 km for the horizontal winds. On the other hand, the MNH-DE1D run has an enhanced variance of the vertical wind but

overestimates the horizontal wind energy between 70 and 10 km. This result has to be appreciated with regard to the fact that MNH-DE1D produces no TKE during the entire simulation, and this is discussed on Section 4. Conversely, the peak in the MNH-DISS run is completely removed, and the split from the high-resolution experiments' horizontal wind spectra occurs at larger scales than MNH-CTRL, at around 40 km. But the discrepancy is the largest with the operational suite AROME, which does not reproduce the slope change at around $\lambda = 60$ km, and an energy surplus between $\lambda = 40$ and 10 km captured by all Meso-NH runs. The AROME spectra also exhibit a truncation at $\lambda = 3\Delta x$, with a drop in variance related to the treatment of orography mentioned in Section 2.3.1, although a peak of energy is still noticeable in the $\lambda \in [4 \text{ km}, 5 \text{ km}]$ range. This peak is nonetheless much weaker than for MNH-CTRL. At the smallest scales of the spectra for ARO-L90 and ARO-L105, a signal on $E_z$ associated with the wavelength $\lambda = 2\Delta x$ seems to indicate the presence of numerical instabilities, but still negligible compared with the peak at 4.5 km. This is probably induced by the high Courant-Friedrichs-Lewy numbers due to the high time step with respect to the high winds at these upper levels of the model. Further investigations are still needed to explain this energetic deficit for the AROME runs in the $\lambda \le 60$ km scales. A possible explanation could be the different numerics used in the AROME model, which is tuned to be more diffusive at upper levels of the model in order to avoid numerical blow-ups that could be critical in operational forecasts. Another important difference concerns the time step: AROME used an operational value of $\Delta t = 50$ s, whereas Meso-NH runs used $\Delta t = 2$ s. Several tests on the dynamics of the AROME model concerning time-step values and the tuning of the numerical solver showed these had a substantial impact on the AROME spectra, but none of the changes that were tested allowed the model to retrieve the same behaviour as the Meso-NH model.

## 3.3 | Characterization of turbulence at small scales

This section considers the properties of the turbulent structures seen on the high-resolution runs.

First, the features of the energy peak in horizontal spectra are discussed, regarding the characteristics of the growing perturbation resolved at $\Delta x = 260$ m. Overturning, occurring between $z = 11$ and 12 km, is signaled with the presence of $Ri_g \le 0$ in the high-resolution simulations. The horizontal wavelengths at which the peak energy exists in the spectra, at $\lambda = 4.5$ km, rather indicate that the scales of the perturbation are small compared with that
<figure><img src="image_20.png" /><figcaption>FIGURE 15 Power density of kinetic energy spectra on domain 3 at 1300 UTC for (a) horizontal wind and (b) vertical wind, averaged from the 175 hPa, 200 hPa, and 225 hPa levels. $κ$: horizontal wave number. [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

<figure><img src="image_21.png" /></figure>

of the mesoscale waves in the ~70 km scale region. As highlighted in Section 2.2, the role of larger scale waves with sufficient amplitude would be to locally modify the properties of thermal stability (Ellrod et al., 2003; Luce et al., 2008; Plougonven and Zhang, 2014), which could then trigger local instabilities. Previous studies also linked perturbations of this scale in the upper flank of a jet to Kelvin-Helmholtz-type instabilities around a critical level (Kim and Chun, 2010; Trier et al., 2020). In the following, we will focus on the characteristics of these instabilities based on linear stability theory.

Wave number and phase velocity of a growing perturbation are determined in the general case through the use of the Taylor-Goldstein equation, relating them to wind and stability profiles of the flow. Since the general solution of this equation is not tractable, normal modes are retrieved using simplified assumptions regarding these profiles (Godrèche and Manneville, 1998; Drazin and Reid, 2004, p. 171; p. 325). Here, we make use of an oversimplified version, assuming $N^2 = 0$, and see how it compares with the model results. In the thermally neutral limit, the wavelength of the fastest growing unstable mode moving at the critical level (i.e., the stationary wave relatively to the mean flow that develops within the shear layer) can be written for piecewise continuous profiles (Cushman-Roisin and Beckers, 2011, p. 327):

$$ \lambda_{\max} = 7.89\delta, \quad (12) $$

where $δ$ is the thickness of the shear layer. Using $δ \sim 500$ m (a value chosen from the layers of $Ri \le 0.25$ in oblique sections, Figure 12, in which the growth of shear instability is expected) gives in our case $λ_{max} = 3.945$ km, a wavelength of the same order of magnitude as the wind disturbances in MNH-HR3D and MNH-HR1D.

The corresponding growth rate from this analysis is $ω_{max} = 0.201ΔU/δ$, with $ΔU$ being the horizontal wind difference across the shear layer. Taking $ΔU = 10 m·s^{-1}$ on the layer considered here gives $ω_{max} = 4.02 × 10^{-3} rad·s^{-1}$, whereas the fit used in Section 3.2.2 gives $1.09 × 10^{-4} ± 7.061 × 10^{-6} rad·s^{-1}$. This is an order of magnitude of difference from the one prescribed from the linear expression in the neutral case, which does not take into account the expense of kinetic energy that the perturbations need to grow against the buoyancy forces (Miles, 1961; Lin, 2007). However, this comparison strongly suggests that the energy peaks represented on the kinetic energy spectra are linked to dynamical, shear instabilities, sub-resolved at $Δx = 1.3$ km.

Since an inertial range scaling seems to be resolved for the spectra of the wind, it should be of interest to further characterize the scales of that turbulence based on classical length scales for stably stratified flows (Lindborg, 2006; Riley and Lindborg, 2008). In particular, the Ozmidov length scale $l_O$ and an outer length scale $l_B$ indicate the range of scales on which we expect to exhibit a kinetic energy cascade. The Ozmidov length indicates the
smallest vertical scale on which buoyancy is expected to impact the turbulent fluctuations and is given by

$$l_O = \left( \frac{\epsilon}{N^3} \right)^{1/2} . \qquad (13)$$

The outer length scale $l_B$ is computed here as an aspect ratio between a characteristic horizontal wind scale $U_h$ and the buoyancy frequency $N$:

$$l_B = \frac{U_h}{N}. \qquad (14)$$

Using values for $z = 11.5$ km, which corresponds to the maximum subgrid TKE layer for MNH-HR3D and MNH-CTRL around 1330 UTC and 1200 UTC respectively, $l_B$ and $l_O$ can be evaluated. Using MNH-CTRL mean values for the large-scale quantities $U_h$ and $N^2$ (Figure 10), which are assumed to vary little between the 2 hr, namely

$$U_h \approx 40 \text{ m} \cdot \text{s}^{-1} \quad \text{and} \quad N^2 \approx 2 \times 10^{-4} \text{ s}^{-2}, \qquad (15)$$

and taking a reference value for $\epsilon$ using the parametrized budget term (Figure 13) from MNH-HR3D between 1300 and 1400 UTC,

$$\begin{gathered} \epsilon \approx 550 \text{ m}^2 \cdot \text{s}^{-2} \cdot \text{day}^{-1} = 6.37 \times 10^{-3} \text{ m}^2 \cdot \text{s}^{-3} \\ \Rightarrow \text{EDR} \approx 1.85 \times 10^{-1} \text{ m}^{2/3} \cdot \text{s}^{-1}, \qquad (16) \end{gathered}$$

gives the following for $l_B$ and $l_O$:

$$l_B \approx 2.82 \text{ km} \quad \text{and} \quad l_O \approx 47.5 \text{ m}. \qquad (17)$$

Thus, according to this estimation, these scales are sufficiently separated and an inertial range can develop. Estimates for $l_B$ and $l_O$ on the same order of magnitudes for CAT have been derived from high-frequency measurements in the study of Rodriguez Imazio et al. (2022). Nonetheless, caution must be emphasized, as it is not clear if the simulations at $\Delta x = 260$ m can be used as a reference. Here, the simulation roughly attains 80% of resolution for the variance on vertical velocity, which is barely at the threshold for LES mode (Pope, 2004).

## 4 | DISCUSSION AND PERSPECTIVES

A case of clear-air jet stream turbulence, selected from available EDR measurements indicating moderate levels of CAT, has been reproduced with the meteorological models AROME and Meso-NH. The large-scale analysis seems to show that this event is mainly forced by the propagation of mesoscale waves in a zone of reduced dynamic

stability, thus triggering the amplification of shear instabilities.

The performance of the models running at operational resolution to reproduce turbulence has been assessed by comparing them with in-situ dissipation measurements. It is shown that enhanced vertical discretization at the tropopause has significant impacts on the EDR forecast, as it allows the better representation of vertical wind shear, which is the leading mechanism that produces subgrid turbulence in the parametrization. However, measurements only give a partial view of the event, given the reduced number of available reports.

The sensitivity of the turbulence to the parametrization and different formulations of the turbulence scheme has been studied. Most notably, the accounting of horizontal gradients in the TKE equation did not enhance subgrid mixing at all during the entire run. This contrasts with the behaviour of turbulence near the ground, where three-dimensional effects are known to have significant impacts on the computed subgrid TKE at the NWP resolution. Also, changing the formulation of the mixing length in the scheme had a strong impact on the ability of the model to produce subgrid TKE in this case. The non-local mixing-length formulations BL89 and RM17 (Bougeault and Lacarrere, 1989; Rodier et al., 2017) succeed at producing subgrid TKE at $\Delta x = 1.3$ km. Finally, the increase of subgrid TKE using different parametrizations also reduces wind and potential temperature gradients, highlighting the diffusive role of the turbulence scheme on the resolved meteorological fields.

High-resolution simulations have also been performed with the Meso-NH mesoscale research model at $\Delta x = 260$ m. They allow one to link the turbulence at the UTLS levels with wind shear instabilities, with the use of wind spectra indicating the characteristic horizontal wavelengths at which they develop. The subgrid budget terms and the subgrid TKE values show a large sensitivity to horizontal gradients at this hectometre resolution. Negative Richardson values inside the shear layer indicate that wave overturning or breaking takes place at the UTLS, which is not represented at the operational resolution. Coarse-graining tools are then used to compare the simulations at $\Delta x = 1.3$ km and 260 m, showing an underestimation of TKE by the current turbulence scheme at $\Delta x = 1.3$ km with the MNH-CTRL run. Also in the wind spectra, the comparison of the different Meso-NH experiments at different resolutions seems to imply that shear instabilities in the jet exit region are only partially resolved by the model at operational resolution with the MNH-DE1D and MNH-CTRL runs, raising further concerns about whether these events are on the so-called “grey zone” (Wyngaard, 2004) of the atmospheric models at $\Delta x = 1.3$ km, where assumptions used in the formulation of
turbulence closures, such as horizontal homogeneity, might have to be revisited (Muñoz-Esparza et al., 2020). Spectra of the MNH-DE1D, MNH-CTRL, and MNH-DISS runs show that runs with increased subgrid mixing in the model also resorb the small-scale fluctuations of the vertical wind at 4.5 km. This points to some transfer of the representation of turbulence between the resolved part of the model and the subgrid values of the scheme in this case. Furthermore, these modified mixing intensities also have a strong impact on the spectra values shown for $\lambda \in [10 \text{ km}, 100 \text{ km}]$, where horizontal wind energy is either overestimated (MNH-DE1D) or underestimated (MNH-DISS) when compared with the high-resolution simulations. Power density spectra were also examined at subsequent hours, confirming the tendencies of these two runs to diverge from the $\Delta x = 260 \text{ m}$ simulations, suggesting that mixing is unphysical (either too strong or too weak) in both cases. A better parametrization should, therefore, be able to satisfy a trade-off between an increase of TKE and other fields from the scheme, which seems to be underestimated in our case, and the alleviation of this induced damping of larger scale fields from the turbulence parametrization. A too large damping could have undesired effects on the representation of other processes at larger scales, such as the depth of low-pressure systems and anticyclones (Sandu et al., 2013).

These questions should be important for modellers, who would have to choose between a jump in resolution to enhance the resolved fields, with its incurring computational cost, or to improve the parametrization of turbulence, in order to enhance the turbulence outputs used as metrics for the operational diagnostics of CAT events. In the latter, high-resolution datasets such as the one presented in this study could be used to evaluate the parametrized quantities of the turbulence scheme at the operational kilometre resolution for the UTLS. Finally, the partial resolution of an inertial range in the kinetic energy spectra would seem to indicate that LES frameworks could be used for the characterization of CAT in realistic set-ups. Further studies should focus on how much spatial resolution is needed for such a task and use these data with the aim to improve the current turbulent parametrizations.

## AUTHOR CONTRIBUTIONS

**Léo Rogel:** data curation; formal analysis; investigation; methodology; software; validation; visualization; writing – original draft. **Didier Ricard:** conceptualization; formal analysis; investigation; methodology; software; super-vision; validation; writing – review and editing. **Éric Bazile:** conceptualization; formal analysis; investigation;

methodology; software; supervision; validation; writing – review and editing. **Irina Sandu:** conceptualization; formal analysis; methodology; supervision.

## ACKNOWLEDGEMENTS

We would like to thank the three reviewers for their very helpful and interesting comments that helped to improve the manuscript in many aspects. We also want to thank Olivier Jaron and Pierre Crispel, who made available the in-situ dataset used in the study. Olivier Jaron also provided the frontogenesis function diagnostic and insights for its interpretation.

## DATA AVAILABILITY STATEMENT

The in-situ observational dataset from the aircraft that has been used for EDR measurements is available on request from the corresponding author.

## ORCID

Léo Rogel https://orcid.org/0000-0003-2146-0623

Didier Ricard https://orcid.org/0000-0002-7323-6300

Irina Sandu https://orcid.org/0000-0002-1215-3288

## REFERENCES

Bazile, E., Marquet, P., Bouteloup, Y. and Bouyssel, F. (2011) The Turbulent Kinetic Energy (TKE) Scheme in the NWP Models at Meteo France. ECMWF, Shinfield Park, Reading: In Workshop on Diurnal Cycles and the Stable Boundary Layer.

Beare, R.J., Macvean, M.K., Holtslag, A.A.M., Cuxart, J., Esau, I., Golaz, J.-C., Jimenez, M.A., Khairoutdinov, M., Kosovic, B., Lewellen, D., Lund, T.S., Lundquist, J.K., Mccabe, A., Moene, A.F., Noh, Y., Raasch, S. and Sullivan, P. (2006) An Intercomparison of large-Eddy simulations of the stable boundary layer. *Boundary-Layer Meteorology*, 118, 247–272.

Belamari, S. and Pirani, A. (2007) Validation of the optimal heat and momentum fluxes using the ORCA2-LIM global ocean-ice model. *Marine EnviRonment and Security for the European Area-Integrated Project (MERSEA IP), Deliverable D*, 4, 550.

Bluestein, H.B. (1992) *Synoptic-Dynamic Meteorology in Midlatitudes: Observations and Theory of Weather Systems*. New York: Taylor & Francis.

Bougeault, P. and Lacarrere, P. (1989) Parameterization of orography-induced turbulence in a Mesobeta-scale model. *Monthly Weather Review*, 117, 1872–1890.

Bramberger, M., Dörnbrack, A., Wilms, H., Ewald, F. and Sharman, R. (2020) Mountain-wave turbulence encounter of the research aircraft HALO above Iceland. *Journal of Applied Meteorology and Climatology*, 59, 567–588.

Brousseau, P., Seity, Y., Ricard, D. and Léger, J. (2016) Improvement of the forecast of convective activity from the AROME-France system. *Quarterly Journal of the Royal Meteorological Society*, 142, 2231–2243.

Bubnová, R., Hello, G., Bénard, P. and Geleyn, J.-F. (1995) Integration of the fully elastic equations cast in the hydrostatic pressure terrain-following coordinate in the framework of the
ARPEGE/Aladin NWP system. Monthly Weather Review, 123, 515–535.

Canuto, V.M., Cheng, Y., Howard, A.M. and Esau, I.N. (2008) Stably stratified flows: a model with No Ri(cr). Journal of the Atmospheric Sciences, 65, 2437–2447.

Colella, P. and Woodward, P.R. (1984) The piecewise parabolic method (PPM) for gas-dynamical simulations. Journal of Computational Physics, 54, 174–201.

Courtier, P., Freydier, C., Geleyn, J., Rabier, F. and Rochas, M. (1991) The ARPEGE project at météo-France. Proceedings of the Numerical Methods in Atmospheric Models, 2, 192–231.

Couvreux, F., Bazile, E., Rodier, Q., Maronga, B., Matheou, G., Chinita, M.J., Edwards, J., van Stratum, B.J.H., van Heerwaarden, C.C., Huang, J., Moene, A.F., Cheng, A., Fuka, V., Basu, S., Bou-Zeid, E., Canut, G. and Vignon, E. (2020) Intercomparison of large-Eddy simulations of the Antarctic boundary layer for very stable stratification. Boundary-Layer Meteorology, 176, 369–400.

Crispel, P. and Jaron, O. (2020) Diagnostic de Turbulence en EDR Pour l'aéronautique Pour ARPEGE cy43t2_op3. Météo-France: Tech. rep.

Cushman-Roisin, B. and Beckers, J.-M. (2011) Introduction to Geophysical Fluid Dynamics: Physical and Numerical Aspects. Waltham: Academic Press.

Cuxart, J., Bougeault, P. and Redelsperger, J.-L. (2000) A turbulence scheme allowing for mesoscale and large-eddy simulations. Quarterly Journal of the Royal Meteorological Society, 126, 1–30.

Deardorff, J.W. (1980) Stratocumulus-capped mixed layers derived from a three-dimensional model. Boundary-Layer Meteorology, 18, 495–527.

Denis, B., Côté, J. and Laprise, R. (2002) Spectral decomposition of two-dimensional atmospheric fields on limited-area domains using the discrete cosine transform (DCT). Monthly Weather Review, 130, 1812–1829.

Dörnbrack, A., Bechtold, P. and Schumann, U. (2022) High-resolution aircraft observations of turbulence and waves in the free atmosphere and comparison with global model predictions. Journal of Geophysical Research: Atmospheres, 127, e2022JD036654.

Dorrestijn, J., Crommelin, D.T., Siebesma, A.P. and Jonker, H.J.J. (2013) Stochastic parameterization of shallow cumulus convection estimated from high-resolution model data. Theoretical and Computational Fluid Dynamics, 27, 133–148.

Drazin, P.G. and Reid, W.H. (2004) Hydrodynamic Stability. Cambridge Mathematical Library, 2nd edition. Cambridge: Cambridge University Press.

Durran, D.R. (1989) Improving the Anelastic approximation. Journal of the Atmospheric Sciences, 46, 1453–1461.

Dutton, J.A. and Panofsky, H.A. (1970) Clear air turbulence: a mystery may be unfolding. Science, 167, 937–944.

Ellrod, G., Lester, P. and Ehernberger, L. (2003) Clear air turbulence. Encyclopedia of Atmospheric Sciences, 3, 393–403.

Fouquart, Y. and Bonnel, B. (1980) Computations of solar heating of the earth’s atmosphere-a new parameterization. Beitraege zur Physik der Atmosphaere, 53, 35–62.

Gal-Chen, T. and Somerville, R.C.J. (1975) On the use of a coordinate transformation for the solution of the Navier-stokes equations. Journal of Computational Physics, 17, 209–228.

Gidel, L.T. and Shapiro, M.A. (1979) The role of clear air turbulence in the production of potential vorticity in the vicinity of upper tropospheric jet stream-frontal systems. Journal of the Atmospheric Sciences, 36, 2125–2138.

Godrèche, C. and Manneville, P. (Eds.). (1998) Hydrodynamics and Nonlinear Instabilities. Collection Alea-Saclay: Monographs and Texts in Statistical Physics. Cambridge: Cambridge University Press.

Goecke, T. and Machulskaya, E. (2021) Aviation turbulence forecasting at DWD with ICON: methodology, case studies, and verification. Monthly Weather Review, 149, 2115–2130.

Goger, B., Rotach, M.W., Gohm, A., Fuhrer, O., Stiperski, I. and Holtslag, A.A.M. (2018) The impact of three-dimensional effects on the simulation of turbulence kinetic energy in a major Alpine Valley. Boundary-Layer Meteorology, 168, 1–27.

Holtslag, A., Svensson, G., Baas, P., Basu, S., Beare, B., Beljaars, A.C.M., Bosveld, F.C., Cuxart, J., Lindvall, J., Steeneveld, G.J., Tjernström, M. and Van De Wiel, B.J.H. (2013) Stable atmospheric boundary layers and diurnal cycles: challenges for weather and climate models. Bulletin of the American Meteorological Society, 94, 1691–1706.

Honnert, R., Masson, V. and Couvreux, F. (2011) A diagnostic for evaluating the representation of turbulence in atmospheric models at the kilometric scale. Journal of the Atmospheric Sciences, 68, 3112–3131.

Hoskins, B.J. (1974) The role of potential vorticity in symmetric stability and instability. Quarterly Journal of the Royal Meteorological Society, 100, 480–482.

Hunt, J.C.R. (1988) Length Scales in Stably Stratified Turbulent Flows and their Use in Turbulence Models, Stably stratified flow and dense gas dispersion.

Kennedy, P.J. and Shapiro, M.A. (1975) The energy budget in a clear air turbulence zone as observed by aircraft. Monthly Weather Review, 103, 650–654.

Kennedy, P.J. and Shapiro, M.A. (1980) Further encounters with clear air turbulence In research aircraft. Journal of the Atmospheric Sciences, 37, 986–993.

Kim, J.-H. and Chun, H.-Y. (2010) A numerical study of clear-air turbulence (CAT) encounters over South Korea on 2 April 2007. Journal of Applied Meteorology and Climatology, 49, 2381–2403.

Kim, J.-H., Chun, H.-Y., Sharman, R.D. and Trier, S.B. (2014) The role of vertical shear on aviation turbulence within cirrus bands of a simulated Western Pacific cyclone. Monthly Weather Review, 142, 2794–2813.

Knox, J.A. (1997) Possible mechanisms of clear-air turbulence in strongly anticyclonic flows. Monthly Weather Review, 125, 1251–1259.

Koch, S.E., Jamison, B.D., Lu, C., Smith, T.L., Tollerud, E.I., Girz, C., Wang, N., Lane, T.P., Shapiro, M.A., Parrish, D.D. and Cooper, O.R. (2005) Turbulence and gravity waves within an upper-level front. Journal of the Atmospheric Sciences, 62, 3885–3908.

Kolmogorov, A.N. (1941) The local structure of turbulence in incompressible viscous fluid for very large Reynolds numbers. The Proceedings of the URSS Academy of Sciences, 30, 301–305.

Lac, C., Chaboureau, J.-P., Masson, V., Pinty, J.-P., Tulet, P., Escobar, J., Leriche, M., Barthe, C., Auizerats, B., Augros, C., Aumond, P., Auguste, F., Bechtold, P., Berthet, S., Bielli, S., Bosseur, F., Caumont, O., Cohard, J.-M., Colin, J., Couvreux, F., Cuxart, J., Delautier, G., Dauhut, T., Ducrocq, V., Filippi, J.-B., Gazen, D., Geoffroy, O., Gheusi, F., Honnert, R., Lafore, J.-P., Lebeaupin Brossier, C., Libois, Q., Lunet, T., Mari, C., Maric, T., Mascart, P., Mogé, M., Molinié, G., Nuissier, O., Pantillon, F., Peyrillé, P.,
Pergaud, J., Perraud, E., Pianezze, J., Redelsperger, J.-L., Ricard, D., Richard, E., Riette, S., Rodier, Q., Schoetter, R., Seyfried, L., Stein, J., Suhre, K., Taufour, M., Thouron, O., Turner, S., Verrelle, A., Vié, B., Visentin, F., Vionnet, V. and Wautelet, P. (2018) Overview of the Meso-NH model version 5.4 and its applications. *Geoscientific Model Development*, 11, 1929–1969.

Laprise, R. (1992) The Euler equations of motion with hydrostatic pressure as an independent variable. *Monthly Weather Review*, 120, 197–207.

Lin, Y.-L. (2007) *Mesoscale Dynamics*. Cambridge: Cambridge University Press.

Lindborg, E. (2006) The energy cascade in a strongly stratified fluid. *Journal of Fluid Mechanics*, 550, 207–242.

Luce, H., Hassenpflug, G., Yamamoto, M., Fukao, S. and Sato, K. (2008) High-resolution observations with MU radar of a KH instability triggered by an inertia–gravity wave in the upper part of a jet stream. *Journal of the Atmospheric Sciences*, 65, 1711–1718.

Mahalov, A., Moustaoui, M., Nicolaenko, B. and Tse, K.L. (2007) Computational studies of inertia-gravity waves radiated from upper tropospheric jets. *Theoretical and Computational Fluid Dynamics*, 21, 399–422.

Mahrt, L. (2014) Stably stratified atmospheric boundary layers. *Annual Review of Fluid Mechanics*, 46, 23–45.

Masciadri, E., Lascaux, F., Turchi, A. and Fini, L. (2017) Optical turbulence forecast: ready for an operational application. *Monthly Notices of the Royal Astronomical Society*, 466, 520–539.

Masciadri, E., Vernin, J. and Bougeault, P. (1999) 3D mapping of optical turbulence using an atmospheric numerical model-I. a useful tool for the ground-based astronomy. *Astronomy and Astrophysics Supplement Series*, 137, 185–202.

Masson, V. (2000) A physically-based scheme for the urban energy budget In atmospheric models. *Boundary-Layer Meteorology*, 94, 357–397.

Meischner, P., Baumann, R., Höller, H. and Jank, T. (2001) Eddy dissipation rates in thunderstorms estimated by Doppler radar in relation to aircraft In situ measurements. *Journal of Atmospheric and Oceanic Technology*, 18, 1609–1627.

Miles, J.W. (1961) On the stability of heterogeneous shear flows. *Journal of Fluid Mechanics*, 10, 496–508.

Miles, J.W. (1963) On the stability of heterogeneous shear flows. Part 2. *Journal of Fluid Mechanics*, 16, 209–227.

Mlawer, E.J., Taubman, S.J., Brown, P.D., Iacono, M.J. and Clough, S.A. (1997) Radiative transfer for inhomogeneous atmospheres: RRTM, a validated correlated-k model for the longwave. *Journal of Geophysical Research*, 102, 16663–16682.

Muñoz-Esparza, D., Sharman, R.D. and Lundquist, J.K. (2018) Turbulence dissipation rate in the atmospheric boundary layer: observations and WRF mesoscale Modeling during the XPIA field campaign. *Monthly Weather Review*, 146, 351–371.

Muñoz-Esparza, D., Sharman, R.D. and Trier, S.B. (2020) On the consequences of PBL scheme diffusion on UTLS wave and turbulence representation in high-resolution NWP models. *Monthly Weather Review*, 148, 4247–4265.

Noilhan, J. and Planton, S. (1989) A simple parameterization of land surface processes for meteorological models. *Monthly Weather Review*, 117, 536–549.

Pergaud, J., Masson, V., Malardel, S. and Couvreux, F. (2009) A parameterization of dry thermals and shallow cumuli for mesoscale numerical weather prediction. *Boundary-Layer Meteorology*, 132, 83.

Pinty, J.-P. and Jabouille, P. (1998) A mixed-phase cloud parameterization for use in mesoscale non-hydrostatic model: simulations of a squall line and of orographic precipitations. In *Conf. On Cloud Physics*, 217–220. Amer. Meteor. Soc Everett, WA.

Plougonven, R. and Zhang, F. (2014) Internal gravity waves from atmospheric jets and fronts. *Reviews of Geophysics*, 52, 33–76.

Pope, S.B. (2004) Ten questions concerning the large-eddy simulation of turbulent flows. *New Journal of Physics*, 6, 35.

Prosser, M.C., Williams, P.D., Marlton, G.J. and Harrison, R.G. (2023) Evidence for large increases in clear-air turbulence over the past four decades. *Geophysical Research Letters*, 50, e2023GL103814.

Redelsperger, J.L. and Sommeria, G. (1981) Methode de representation de la turbulence d'echelle inferieure a la maille pour un mod- ele tri-dimensionnel de convection nuageuse. *Boundary-Layer Meteorology*, 21, 509–530.

Reiter, E.R. (1969) The nature of clear air turbulence: a review. In: Pao, Y.-H. and Goldburg, A. (Eds.) *Clear Air Turbulence and its Detection*. Boston, MA: Springer US, pp. 7–33.

Ricard, D., Lac, C., Riette, S., Legrand, R. and Mary, A. (2013) Kinetic energy spectra characteristics of two convection-permitting limited-area models AROME and Meso-NH. *Quarterly Journal of the Royal Meteorological Society*, 139, 1327–1341.

Riley, J.J. and Lindborg, E. (2008) Stratified turbulence: a possible interpretation of some geophysical turbulence measurements. *Journal of the Atmospheric Sciences*, 65, 2416–2424.

Rodier, Q., Masson, V., Couvreux, F. and Paci, A. (2017) Evaluation of a buoyancy and shear based mixing length for a turbulence scheme. *Frontiers Earth Science*, 5.

Rodriguez Imazio, P., Dörnbrack, A., Urzua, R.D., Rivaben, N. and Godoy, A. (2022) Clear air turbulence observed across a tropopause fold over the Drake Passage–a case study. *Journal of Geophysical Research: Atmospheres*, 127, e2021JD035908.

Sandu, I., Beljaars, A., Bechtold, P., Mauritsen, T. and Balsamo, G. (2013) Why is it so difficult to represent stably stratified conditions in numerical weather prediction (NWP) models? *Journal of Advances in Modeling Earth Systems*, 5, 117–133.

Seity, Y., Brousseau, P., Malardel, S., Hello, G., Bénard, P., Boutier, F., Lac, C. and Masson, V. (2011) The AROME-France convective-scale operational model. *Monthly Weather Review*, 139, 976–991.

Shapiro, M.A. (1976) The role of turbulent heat flux in the generation of potential vorticity in the vicinity of upper-level jet stream systems. *Monthly Weather Review*, 104, 892–906.

Sharman, R. and Lane, T. (Eds.). (2016) *Aviation Turbulence: Processes, Detection*. Prediction: Springer International Publishing.

Sharman, R., Tebaldi, C., Wiener, G. and Wolff, J. (2006) An integrated approach to mid- and upper-level turbulence forecasting. *Weather and Forecasting*, 21, 268–287.

Sharman, R.D., Cornman, L.B., Meymaris, G., Pearson, J. and Far-rar, T. (2014) Description and derived Climatologies of automated In situ Eddy-dissipation-rate reports of atmospheric turbulence. *Journal of Applied Meteorology and Climatology*, 53, 1416–1432.

Sharman, R.D., Doyle, J.D. and Shapiro, M.A. (2012a) An investigation of a commercial aircraft encounter with severe clear-air turbulence over Western Greenland. *Journal of Applied Meteorology and Climatology*, 51, 42–53.

Sharman, R.D. and Pearson, J.M. (2016) Prediction of energy dissipation rates for aviation turbulence. Part I: forecasting nonconvective turbulence. *Journal of Applied Meteorology and Climatology*, 56, 317–337.
<figure><img src="image_22.png" /></figure>

Sharman, R.D., Trier, S.B., Lane, T.P. and Doyle, J.D. (2012b) Sources and dynamics of turbulence in the upper troposphere and lower stratosphere: a review. *Geophysical Research Letters*, 39.

Skamarock, W., Klemp, J., Dudhia, J., Gill, D., Barker, D., Wang, W., Huang, X.-Y. and Duda, M. (2008) *A Description of the Advanced Research WRF Version 3*. UCAR/NCAR: Tech. rep.

Skamarock, W.C., Park, S.-H., Klemp, J.B. and Snyder, C. (2014) Atmospheric kinetic energy spectra from global high-resolution nonhydrostatic simulations. *Journal of the Atmospheric Sciences*, 71, 4369–4381.

Strauss, C., Ricard, D. and Lac, C. (2022) Dynamics of the cloud-Environment Interface and turbulence effects in an LES of a growing cumulus Congestus. *Journal of the Atmospheric Sciences*, 79, 593–619.

Trier, S.B., Sharman, R.D., Muñoz-Esparza, D. and Lane, T.P. (2020) Environment and mechanisms of severe turbulence in a Midlatitude cyclone. *Journal of Atmospheric Sciences*, 77, 3869–3889.

Uccellini, L.W. and Koch, S.E. (1987) The synoptic setting and possible energy sources for mesoscale wave disturbances. *Monthly Weather Review*, 115, 721–729.

Verrelle, A., Ricard, D. and Lac, C. (2015) Sensitivity of high-resolution idealized simulations of thunderstorms to horizontal resolution and turbulence parametrization. *Quarterly Journal of the Royal Meteorological Society*, 141, 433–448.

Wacławczyk, M., Gozingan, A.S., Nzotungishaka, J., Mohammadi, M. and Malinowski, S.P. (2020) Comparison of different techniques to calculate properties of atmospheric turbulence from low-resolution data. *Atmosphere*, 11, 199.

Williams, P.D. (2017) Increased light, moderate, and severe clear-air turbulence in response to climate change. *Advances in Atmospheric Sciences*, 34, 576–586.

Wyngaard, J.C. (2004) Toward numerical Modeling in the “Terra incognita”. *Journal of the Atmospheric Sciences*, 61, 1816–1826.

Zhang, F., Koch, S.E., Davis, C.A. and Kaplan, M.L. (2000) A survey of unbalanced flow diagnostics and their application. *Advances in Atmospheric Sciences*, 17, 165–183.

**How to cite this article:** Rogel, L., Ricard, D., Bazile, É. & Sandu, I. (2023) Effects of subgrid-scale turbulence parametrization on the representation of clear-air turbulence using kilometre- to hectometre-scale numerical simulations. *Quarterly Journal of the Royal Meteorological Society*, 149(757), 3301–3322. Available from: https://doi.org/10.1002/ qj.4557