A Theoretical Analysis of Mixing Length for Atmospheric Models From Micro to Large Scales

A new mixing length adapted to the constraints of the hectometric-scale gray zone of turbulence for neutral and convective boundary layers is proposed. It combines a mixing length for mesoscale simulations, where the turbulence is fully subgrid and a mixing length for Large-Eddy Simulations, where the coarsest turbulent eddies are explicitly resolved. The mixing length is built for isotropic turbulence schemes, as well as schemes using the horizontal homogeneity assumption. This mixing length is tested over three boundary layer cases: a free convective case, a neutral case and a cold air outbreak case. The later combines turbulence from thermal and dynamical origins as well as presence of clouds. With this new mixing length, the turbulence scheme produces the right proportion between subgrid and resolved turbulent exchanges in Large Eddy Simulations, in the gray zone and at the mesoscale. This opens the way of using a single mixing length whatever the grid mesh of the atmospheric model, the evolution stage or the depth of the boundary layer.


INTRODUCTION
The atmospheric boundary layer (ABL) is the part of the atmosphere directly in contact with and rapidly influenced by the surface (Stull, 1988). This region is characterized by its turbulence. Wyngaard (2004) defines the gray zone of turbulence as the grid scale close to the scale of the energy containing structures. These structures characteristics depend on the ABL type, convective or neutral, for instance. Honnert et al. (2011) studied the resolved and subgrid parts of dry and cumulus-topped convective ABL (denoted hereafter CBL) simulations between 62.5 m grid spacing (i.e., large-eddy simulations or LES) and 8 km grid spacing (i.e., mesoscale regime). They found that the gray zone of turbulence ranges between 0.2 and 2 times the ABL height in convective cases (or the top of the cloud layer in cumulus-topped cases). Beare (2014) has defined an effective length scale for numerical models, l d,eff , accounting for the modeled energy dissipation emerging from the subgrid and the advection scheme. He showed that the transition between the gray zone and the mesoscale occurs at h/l d,eff 0.7 in a convective case (where h is the boundary-layer height). Honnert (2018), using the method of Honnert et al. (2011), showed that the gray zone ranges between 0.03h and h in neutral cases. In a nutshell, the gray zone of turbulence concerns the hectometric scales.
Numerical Weather Prediction (NWP) models are nowadays operationally used at the mesoscale (kilometric scale) and need improved turbulence schemes in the gray zone. At the mesoscale, large scale turbulence may be represented by non-local diffusion, using a non-local Eddy-diffusivity-based scheme (Lock et al., 2000) or a mass-flux scheme (Hourdin et al., 2002;Siebesma et al., 2004;Pergaud et al., 2009). However, Honnert et al. (2011) showed that these turbulence schemes may have difficulties to adapt in the gray zone of turbulence. Typical mesoscale turbulence schemes such as mass-flux schemes (e.g., PM09, Pergaud et al., 2009) produce too much subgrid mixing. Honnert et al. (2016) and Lancz et al. (2017) proposed to weaken the mass-flux scheme in order to produce resolved turbulent structures in the gray zone of turbulence. However, at the hectometric scales, non-local structures may be widely resolved whereas the rest of the turbulence, the local turbulence, cannot be considered as homogeneous and isotropic and thus cannot be directly treated by current LES schemes (Honnert and Masson, 2014). As a consequence, the local turbulence has also to be adapted to the gray zone.
The local turbulence in meteorological models is classically represented using eddy-diffusivity turbulence schemes. In a Smagorinsky scheme (Smagorinsky, 1963), for instance, the eddy-diffusivity depends on the shear and the mixing length.The latter usually represents the typical size of the energy-containing subgrid eddies. Boutle et al. (2014) were the first to build a seamless turbulence scheme in the gray zone. They pragmatically blended a non-local scheme (Lock et al., 2000) and a 3D Smagorinsky scheme (Smagorinsky (1963), appropriate for well-resolved large eddies in LES regime) in a way that is scale-aware and related to the ratio of the grid-scale to a diagnosed length scale for the turbulence. The parameterization requires to modify the computation of the mixing length but also the subgrid flux blending the fluxes from non-local and 3D Smagorinsky schemes. This blending is effective even when it is not necessary: when the non-local turbulence is entirely resolved, fully 3D, or mainly horizontal, which happens in the gray zone at fine resolutions or over complex terrains.
Similarly, Ito et al. (2015) proposed an extension of Mellor-Yamada-Nakanishi-Niino (MYNN) model where the length scale is modified in order to maintain the turbulence kinetic energy (TKE) dissipation invariant to the grid spacing. The modification uses the partition function proposed by Honnert et al. (2011) in order to provide an appropriate subgrid TKE contribution. By using a coarse-graining approach on LES data, Kitamura (2015) estimated the length scale dependency on the grid spacing at the hectometric scale. He assumed that the turbulent fluxes have the form of a TKE-based model (Deardorff, 1980) and proposed two scale-aware mixing lengths, one vertical and the other horizontal. Recently, Kurowski and Teixeira (2018) defined a new mixing length by merging the existing mixing lengths formulations from LES and NWP model. This mixing length may 'overshoot': when the large scale and the LES mixing lengths are close one to the other (which may happen in the gray zone), the resulting mixing length may be larger than each of the two mixing lengths.
At the same time, Efstathiou et al. (2018) worked on the dissipation of the Smagorinsky scheme at fine resolutions by adapting the mixing length to the growth of a CBL. The turbulence spin-up was improved in simulations where the resolution was coarser than in LES. However, the technique seems too expensive for operational applications at the hectometric scales. For the coarser resolutions of the gray zone of turbulence,Ďurán et al. (2020) modified the Bougeault and Lacarrère (1989) mixing length based on budget of TKE for coarsegrained LES. However, they did not intend to reach the LES scales.
NWP models cover a range of scale from 300 m to the mesoscale (2 km and above). They could even reach 100 m resolution in the foreseeable future. For this range of scale, turbulent eddies in the stable boundary layers are still entirely subgrid and thus will not be treated in this article. Indeed, even for a weakly stable boundary layer, a 3vm grid spacing (or even less) is necessary to be considered resolved (Beare et al., 2006). Furthermore, Beare et al. (2006) have shown the sensitivity to Smagorinsky mixing length, backscatter or dynamic model at a grid length of 6.25 m. The sensitivity to subgrid model can be considered as a hallmark of the gray zone. Although, authors dealing with the stable boundary layer gray zone do not call it the gray zone, just an unresolved LES, the gray zone of turbulence in stable boundary layers can be handled by stochastic parameterizations (Brown et al., 1994). or dynamic (Esau, 2004;Khani and Waite, 2014;Khani, 2018). However, depending on the grid-spacing and the boundary layer depth, models do represent eddies in neutral and convective boundary layers either as LES, gray zone or uniquely vertical subgrid motions. Most of the previous-mentioned attempts to modify the turbulence schemes in the gray zone of turbulence are restricted to a part of the scale range from LES to mesoscale.
The objective of this paper is then to propose a new unique mixing length adapted to the gray zone of turbulence for both convective and neutral ABL, pragmatic but based on theoretical foundations and allowing a robust use from LES (typically 50-100 m) to NWP models (at any scale). The paper is organized as follows. A new adaptive formulation for l in the gray zone of turbulence is presented in Section 2. This new mixing length is tested over three boundary layer cases. Their descriptions and the evaluation method are detailed in Section 3. Results are presented in Section 4. A discussion and some future directions are discussed in Section 5. At the end, a summary of findings is given in Section 6.

Mesoscale Mixing Lengths
In climate, NWP or mesoscale research models, non-local mixing lengths are commonly used. Figure 1 represents a schematic eddy created by a non-local mixing length in a large scale model. At large scales, the horizontal resolution (larger than 1 km) is coarser than the vertical resolution closed to the surface (usually finer than a few tens of meters). A non-local mixing length represents large ABL eddies included in the horizontal grid but larger than the vertical one. Non-local mixing lengths have been developed to represent coherent structures, particularly for convective cases (Bougeault and Lacarrere, 1989;Holtslag and Boville, 1993;Lenderink and Holtslag, 2004). For instance, Bougeault and Lacarrère (1989) (hereafter named as BL89) is used in operational NWP models such as AROME (Seity et al., 2011) or WRF 1 (Shin and Hong, 2016).
The size of the largest ABL eddies in BL89 is given by: where l up and l down are the respective mixing lengths for upward and downward movements, z is the altitude, θ v is the virtual potential temperature and β is the buoyancy factor. The mixing length is: Local formulations of the mixing lengths can be based on the local vertical wind shear (Hunt, 1988;Tjernström, 1993;Schumann and Gerz, 1995;Grisogono and Enger, 2004;Grisogono and Belušić, 2008;Grisogono, 2010;Venayagamoorthy and Stretch, 2010). Recently, Rodier et al. (2017) (hereafter RM17) have improved BL89 to make it physically more accurate for all stratification conditions. Indeed, BL89 is infinite for neutral cases. However, in case of strong wind shear, the ABL turbulence is disorganized and the size of the turbulent eddies is reduced compared with a freeconvective case. In order to reproduce the shear-induced vertical disorganization of the free-buoyancy eddies, RM17 integrates the local vertical wind shear within the original BL89 formulation. The size of the largest ABL eddies in RM17 becomes: where C 0 is a constant and σ z||U|| zxj is the local vertical wind shear.
The mixing length is then: Redelsperger et al. (2001)correction is applied in the first 25 m of the model to improve the representation of the wind profile in the surface boundary layer. Moreover, Rodier et al. (2017) promoted the use of a dissipation constant (C ϵ ) reduced from 0.85 (in Bougeault and Lacarrère, 1989) to 0.34.

Large-Eddy Simulation Mixing Lengths
Large-Eddy simulations (Lilly, 1967;Deardorff, 1980;Pope, 2000) resolve most of the turbulence scales. LES are commonly used as a state-of-the-art numerical tool in order to study turbulence or to develop and evaluate parameterizations Rio et al., 2010).
In LES, the grid is much smaller than the largest turbulent eddies and the nearly cubic grid fits the residual isotropic subgrid turbulence (see Figure 2). In this case, the most energycontaining subgrid eddies have the size of the grid box: where Δx, Δy and Δz are the grid cell dimensions. This mixing length (hereafter L DELT ) is prescribed in most LES. However, close to the surface or in very stably stratified cases, the most energy-containing subgrid eddies may be smaller than the grid cell size. This is why Deardorff (1980) mixing length (hereafter DEAR) takes the stratification into consideration: where N is the Brunt-Väïsälä frequency.

A New Mixing Length From Micro to Large Scales
None of the above-described mixing lengths are adapted to the gray zone of turbulence. Indeed, L RM17 or L BL89 hardly adapt to the horizontal resolution, a crucial element in the gray zone when subgrid turbulence mixing must be smaller than the one at large scale. Concerning L DELT and L DEAR , they both depend on the vertical resolution which limits the turbulence mixing, as in Figure 3, whereas subgrid non-local structures exist in the gray zone of turbulence (Shin and Hong, 2014;Honnert et al., 2016). With L DELT , the mixing length formulation takes into account the vertical resolution Δz. In meteorological models, Δz is smaller than the horizontal resolution in Eq. (4) and limits L DELT (see Figure 3). We propose to use a non-local mixing length L , which combines L RM17 with a LES-adapted mixing length L Δ (Eq. 6) that avoids the purely numerical dependency of the mixing length to the model vertical discretization: L is the minimum mixing length between the horizontal grid cell L Δ and L RM17 : where α is a tuning parameter set to 0.5, according to the numerical experiments detailed in Section 4.3. This formulation is continuous and does not overshoot. Usually, L RM17 is of the order of the boundary layer height (but smaller). L is equal to L RM17 for large horizontal grid meshes. At the surface, L RM17 is smaller than L Δ as it is limited by the distance to the surface. For LES (cubic) grid meshes, L is equivalent to Deardorff: the base mixing length is the horizontal grid mesh and it is limited by a stability-based length (RM17), as done in Deardorff length by the TKE part (see Eq. 5) but taking into account shear as well. For grid meshes in the gray zone of turbulence, this is the smaller of the two. The L mixing length only depends on the horizontal grid mesh (which is a key element to define where we stand in regards to the gray zone) and the atmospheric conditions (buoyancy and shear). L is independent of the vertical grid chosen.

Model
The model used in the present work is the atmospheric non hydrostatic research model Meso-NH (Lac et al., 2018). It is used here from LES to hectometric scales in idealized frameworks. The spatial discretization is based on the C grid of Arakawa with Cartesian coordinates. The advection scheme for momentum variables is a centered scheme of 4th order in space, and a Runge-Kutta time-splitting of 4th order in time, well adapted to LES as the effective resolution is 5 − 6Δx (Lac et al., 2018). For potential temperature and hydrometeor mixing ratios, the advection scheme is PPM (Piecewise Parabolic Method) from Colella and Woodward (1984). A numerical diffusion is added to filter the shortest wavelengths and avoid energy accumulation in the tail of the energy spectrum. It is set using a characteristic time (e-folding time) applied to the momentum variables. The atmospheric model is coupled with the EXternalized SURFace (SURFEX) model (Masson et al., 2013). The turbulence scheme is based on CBR (Cuxart et al., 2000), which is a 1.5 order turbulence parameterization based on a TKE prognostic equation, where the eddy-diffusivity, K ϕ , depends on the TKE and a mixing length l:  Frontiers in Earth Science | www.frontiersin.org January 2021 | Volume 8 | Article 582056 where c ϕ is a constant depending of ϕ and e is the TKE. At mesoscale resolutions (horizontal mesh larger than 1 km), it can be assumed that the horizontal gradients and the horizontal turbulent fluxes are much smaller than their vertical counterparts: they are therefore neglected (except for the advection of TKE) and the turbulence scheme is used in its 1D version, as in AROME (Seity et al., 2011). At finer resolution, the whole subgrid equation system is considered in its 3D version, as it is the case in the present work. The microphysical scheme is a mixed-phase one moment scheme (Pinty and Jabouille, 1998) which considers three solid (ice crystal, snow, and graupel) and two liquid (cloud droplets and raindrops) hydrometeors. The radiation parameterizations come from ECMWF, and include the Rapid Radiation Transfer Model (Mlawer et al., 1997) for long-wave and Fouquart and Bonnel (1980) for shortwave radiations.

Case Studies and LES Specifications
Boundary-layer turbulence comes from boundary-layer convection (thermal-induced turbulence) and wind shear (dynamical-induced turbulence). Three cases combining differently one or the other source of turbulence are evaluated here, in dry or cloudy conditions. For all cases, the reference LES is run with L DEAR at high resolution.

Free Convective Boundary Layer
The first case is a free convective boundary layer (CBL). It is based on radio-soundings collected during the International H 2 O project (IHOP 2002 ) campaign (Weckwerth et al., 2004). In this work, we reproduce the June 14th, 2002 case, corresponding to a growing CBL near Homestead, Oklahoma (Couvreux et al., 2005). This day is characterized by high pressure (1,016 hPa or more) and light wind (＜5 m s −1 ). The weak vertical shear guarantees a free convective case. The well-mixed boundary layer reaches 1.2 km in the beginning of the afternoon. This case is chosen as it presents a relatively uniform site topography and a development of a typical continental CBL. The numerical domain is a 16 × 16 km 2 area. Lateral boundary conditions are periodic. Seven hours of dynamics are reproduced. The horizontal grid mesh for the LES is 50 m. The vertical grid is stretched, but the spacing is always ＜100 m in the ABL, so that the grid boxes are close to cubes. The numerical diffusion is small, the time at which the 2Δx waves are damped by the factor e −1 is equal to 1800s. The diffusion value is identical for all resolutions and tested mixing lengths, allowing to compare the results.

Neutral Boundary Layer
The second case is a neutral boundary layer corresponding to the incoming flow conditions of case 2681829 of the Mock Urbain Setting Test experiment (hereafter must, Yee and Biltoft, 2004). This near full-scale experiment was conducted during September 2001 in Utah's West desert, at US Army Dugway Proving Ground. The site is located 1310 m above the mean sea level and can be considered as flat and homogeneous. MUST, starting the September 25, 2001 at 1830 LT (local time UTC−6 h), is characterized by a neutral state atmosphere, a strong wind ( ≈ 8 m s −1 at z 4 m height) and a shear-induced turbulence. This case is chosen as it represents a typical continental neutral boundary layer. Moreover, this case is of particular interest for further developments of the Immersed Boundary Method (IBM) version of the Meso-NH (MNH-IBM, Auguste et al., 2019). The numerical domain is a 40 km side square. The vertical elevation of the domain goes up to 1500 m high. The first vertical layer measures 10 m high. The vertical grid size increases linearly with a geometric ratio of 1.09 and reaches it maximum value (Δz 50 m) at the elevation of 500 m. The vertical grid spacing remains constant between 500 m and the top of the domain. The ground friction is characterized by a roughness length z 0 0.045 m and modeled with SURFEX. An imposed geostrophic wind maintains the flow. Near the ground, the wind conditions match the observations: the wind speed is around 9 m s −1 at z 10 m and the wind direction has an angle of −41°with respect to the x-axis. The case is purely neutral. The L mixing length is tested at 50, 100, 200 and 400 m horizontal resolutions and compared with the reference LES at 50 m. The sensitivity of the results to α is also investigated. The tuning of the numerical diffusion is the same as for the IHOP case.

Cold Air Outbreak
The third case is a cold air outbreak. Cold-air outbreaks are largescale cold air masses over relatively warm waters. This phenomenon is common in the North Atlantic during winter (Field et al., 2014). The CONSTRAIN campaign (Field et al., 2014) studied cold air outbreaks in the North Sea involving a complex microphysics with mixed phases of ice and liquid water and mesoscale 2D (cloud streets) or 3D (open cells) structures. Several models inter-comparisons have been conducted on this study-case: an inter-comparison of global (Tomassini et al., 2016), limited area models (Field et al., 2017) and LES simulations (de Roode et al., 2019). The present case is described in de Roode et al. (2019). It corresponds to a transition from stratocumulus to cumulus as the cloud deck is advected over an increasingly warmer sea surface temperature. This simulation combines turbulence from thermal and dynamical origins, as well as presence of clouds. Similarly to de Roode et al. (2019)

Simulations at Gray Zone Resolutions and Their Evaluation
According to Honnert et al. (2011), the gray zone of turbulence ranges between 0.2 and 2 times the boundary-layer height for CBL. The boundary-layer height being about 1.5 km for both IHOP and CONSTRAIN experiments, horizontal resolutions between 100 and 500 m are therefore considered as in the gray zone of turbulence. For neutral cases, Honnert (2018) has shown that the gray zone of turbulence ranges between 0.03 and 1 times the boundary-layer height. For MUST, the boundary layer is of 1500 m in height and horizontal resolutions between 50 and 400 m are considered as in the gray zone of turbulence.
In kilometric-scale simulations, non-local schemes such as mass-flux schemes (Hourdin et al., 2002;Siebesma et al., 2004;Pergaud et al., 2009) are used to represent the non-local turbulence of CBLs. However, Honnert et al. (2011) proved that mesoscale mass-flux schemes are too diffusive in the gray zone of turbulence and cannot be used at hectometric scales. In addition, the presence of a non-local turbulence scheme would distort the new mixing length impact analysis on the turbulence scheme. This is the case even if the influence of the non-local scheme would be very reduced as in non-local schemes dedicated to the gray zone (Boutle et al., 2014;Honnert, 2016). In the present paper, all simulations are therefore undertaken without any additional nonlocal boundary-layer scheme. We will focus on resolutions between LES and 500 m (upper gray-zone). The L mixing length is assumed to naturally work at mesoscale and large scale because it is then equal to the uni-dimensional mixing length L RM17 . This study does not use direct observations. However, as seen in the previous section, the test cases come from intercomparison studies, meaning that they are close but simplified compared with the reality. We have respected the intercomparison setup where the LES simulations have been validated compared with the observations, allowing the use of those LES as a reference for the evaluation of the new mixing length simulations.
Quantification of the resolved and subgrid fluxes are key points for the parameterization evaluation in the gray zone. Honnert et al. (2011) use the coarse-graining of LES fields in order to quantify the 'true' resolved and subgrid part of the turbulence in the gray zone. It appears that the partitioning between subgrid and resolved fluxes in the mixed layer of CBLs follows one unique function depending only on ΔxΔy √ h . This function is quantified in Honnert et al. (2011) and called partial similarity function. In this work, simulations with L, L DEAR and L RM17 in the gray zone are compared with LES results and partial similarity functions. Figure 4 shows the boundary-layer height computed from the IHOP simulations at resolutions ranging from 50 to 500 m. Three different mixing lengths are used: L, L DEAR and L RM17 . The boundary-layer height is computed as the maximum of the potential temperature gradient during the 7 h simulated.

Boundary-Layer Height
During the first 5 h of dynamics and for all the configurations, the boundary layer height growth rate is similar to the reference LES (dark triangles in Figure 4) but with slightly lower values. Then, the growth of the CBL slows in the LES, which is not the case in high resolution simulations using L RM17 and L DEAR . Table 1 presents the bias and standard deviation of the simulations compared with the LES. Simulations using L have the smallest bias at resolutions smaller than 200 m whereas simulations using L RM17 have the smallest bias at resolutions coarser than 200 m. If the boundary layer develops correctly whatever the configuration, simulations using L are in better agreement with the LES than those using L DEAR . In particular, L DEAR produces a too low boundary-layer height at hectometric scales in agreement with Honnert et al. (2011).

Partial Similarity Functions
At the mesoscale, the turbulence is entirely subgrid and the boundary-layer eddies are smaller than the resolution of the model. In LES, the turbulence is mainly resolved, as it explicitly represents the largest and most energetic turbulent  eddies. In-between, in the gray-zone of turbulence, the subgrid/resolved partitioning of the TKE follows the partial similarity function from Honnert et al. (2011). Figure 5 presents this subgrid/resolved TKE partitioning for the different IHOP simulations undertaken as a function of Δx/h. In the mixed layer of a clear-sky convection boundary-layer, the subgrid/resolved partitioning only depends on the model resolution (whose proxy is the model horizontal grid spacing) and on the size of the largest boundary-layer eddies (whose proxy is the boundary-layer height, assuming that, the higher the boundary-layer eddies, the coarsest they are). More details can be found in Honnert et al. (2011) or Shin and Hong (2014). The black line in Fig. 5 is the partial similarity function of the TKE. The dashed lines are the first and the last vigintiles of the coarse-grained LES data. They represent the allowed margin of error. The data obtained with simulations using L RM17 , L DEAR and L as mixing length are shown in blue, green and red, respectively. L RM17 produces a correct subgrid/resolved partitioning for boundary-layer height smaller than the grid spacing Δx h > 1 . Otherwise Δx h < 1 , the subgrid part of the TKE drops at about 0.6 times the total TKE. This is particularly unrealistic in LES where the TKE should be resolved when Δx h < 0.05. L RM17 presents also inconsistency in the gray zone: the TKE is either too subgrid for Δx h < 0.2, or too resolved for Δx h > 0.3. On the contrary, simulations with L DEAR correctly follow the partial similarity function when Δx h < 0.1 but are always too resolved in the gray zone. Finally, simulations using L follow the partial similarity function until Δx h ≈ 0.4 but are too resolved for 0.4 < Δx h < 0.5. This may be due to the fact that a non-local part of the turbulence is still necessary for resolutions between 400 and 500 m (see Honnert et al., 2016). It has been found (not shown in Figure 5) that the subgrid TKE evolution as a function of Δx h presents very few differences between BL89 and RM17. The impact of the wind shear on turbulence is generally quantified by the −h/L MO ratio, where L MO is the Monin-Obukhov length. For moderate and strong winds, if −h/L MO > 10 the dynamical production is negligible. In IHOP, −h/L MO is close to 100 (see Honnert et al., 2011) for more details) explaining why BL89 and RM17 show similar results.

Spectral Analysis
In order to complete this analysis, Figure 6 presents the energy spectra after 6 h of dynamics, using L or L DEAR . The spectra are computed from data over all the levels between the surface and the top of the boundary layer. Simulations using L RM17 are not shown because whatever the resolution, the TKE is not resolved enough to perform a spectral analysis.
In Fig. 6, the gray dashed line is the Kolmogorov −5/3 slope. The full and dashed lines represent the energy spectra undertaken with L DEAR and L, respectively. There is no significant difference between DEAR and the new mixing length at the finest resolutions. At 50 m grid spacing, L produces energy spectra in agreement with DEAR, considered as the LES reference. The most significant difference is found at 200 m grid spacing where it is clear that DEAR produces too much resolved turbulence compared with the LES and the Kolmogorov theoretical references. L improves the representation of the turbulence at this resolution.

Neutral Boundary Layer: Resolved TKE Vs. Coarse-Grained LES
The neutral MUST case allows to examine the behavior of the mixing lengths in the fine resolution part of the gray zone of  Focusing on the coarse-grained LES results, the TKE is more resolved at the finer resolutions. For instance, TKE is twice more resolved at 100 than at 400 m grid spacing. For all Frontiers in Earth Science | www.frontiersin.org January 2021 | Volume 8 | Article 582056 8 resolutions, no resolved TKE is obtained with RM17. In this neutral case, where sources of instabilities are smaller than in convective cases, the subgrid mixing of RM17 inhibits the resolved movements, even at finest resolutions. This highlights the need of using an appropriate mixing length in the gray zone.
At 50 m grid spacing, the good agreement between LES and L energy profiles for the resolved (Figure 7b) and total TKE (Figure 7a) shows that the new mixing length can be used at high resolution and does not present inconsistency compared with the reference LES. For resolved TKE at coarser resolutions (100 m, Figure 7d, 200 m Figure 7f), L is almost similar, but more resolved than DEAR. The agreement compared with the coarse-grained LES is satisfactory especially at 100 m grid spacing. At 400 m grid spacing (Figure 7h), both L and DEAR overestimate the resolved TKE above 1,000 and 750 m, respectively and underestimate it below. This underestimation is more pronounced with L, especially close to the ground and DEAR is in better agreement with the coarse-grained LES. Nonetheless, L agreement with the coarse-grained LES remains satisfactory in the fine resolution part of the gray zone of turbulence.
The results are similar for the total TKE (Figs 7c, 7e, 7g). L and DEAR are close one to the other and in much better agreement with the LES results than RM17.

Sensitivity Tests to α
The new mixing length L aims to model the turbulent mixing whatever the resolution or the meteorological conditions within the boundary layer. Therefore, the tuning coefficient α should not vary. However, it has to be set. This setting could depend on the atmospheric model in which it is used and more particularly of its intrinsic diffusion properties. In the present work, we will find a compromise to set the value of α for CBR turbulence scheme using the Meso-NH model.
The effective resolution of a model (Skamarock, 2003;Ricard et al., 2013) is the minimum wavelength correctly seen by the model. Due to the physical and numerical diffusion, it is always larger than the grid size (Δx). This effective resolution depends on several parameters: the model considered (more particularly the chosen dynamical scheme), the physical schemes, the internal adjustments of these different schemes and finally the meteorological situation itself. For example, the turbulent diffusion is governed by the chosen mixing length, but also by the other scheme coefficients, such as the turbulence dissipation coefficient. The more diffusive a model, the greater the effective resolution compared to Δx. Meso-NH with the numerical schemes used in this study presents an effective resolution of 5 Δx (Lac et al., 2018) which is the signature of a good accuracy of the numerical schemes. Figure 6 confirms an effective resolution less than 5 Δx, i.e. the scale from which the model departs from the theoretical slope of −5/3. The effective resolution of the model is a main characteristic of the resolved/subgrid partition of the simulations. The theoretical subgrid/resolved partition of Honnert et al. (2011) is obtained from a coarse-graining which gives results from the LES perfect fields. To adjust the parametrization of the mixing length to the diffusion characteristics of the model, a tuning parameter is needed, that is α in Eq. 7.
In IHOP and MUST, the value of α was set to 0.5. The sensitivity of these cases results to α is here investigated. Figure 8 shows the subgrid TKE as a function of Δx h for IHOP (left) and MUST (right) simulations undertaken with α ranging from 0.33 to 1. The black line is the TKE reference partial similarity function from Honnert et al. (2011). The uncertainty margins are represented by the black dashed lines.
Concerning IHOP, Figure 8 left shows that, for the small values of Δx h (i.e., LES regime), large values of α result in diffusive simulations. The subgrid turbulence is too large and the simulation never enters the LES regime because the turbulence is not resolved enough. At the meso-scale (i.e., for large values of Δx h ) the turbulence is entirely subgrid for all values of α. The mixing length is used in the computation of the fluxes, but also in the parameterization of the diffusion. This is particularly obvious in the gray zone, where the larger α, the more diffusive the result is. Moreover, for α 0.33 there is not enough diffusion so that resolved turbulence structures appear even at a 2 km grid spacing ( Δx h ≈ 1.3), which is in contradiction with the reference partial similarity function. Extreme values of α (α 0.33 and α 1) show also a large variability around Δx h 0.8 which is due to waves development in the simulations.
For the MUST case (see Figure 8 right), similarly to IHOP in LES regimes and in the gray-zone, large values of α (α 0.66; 1) result in more diffusive simulations. If subgrid turbulence is small enough to maintain simulations in the LES regimes for Δx h < 0.1, this is not the case for larger Δx h values and results are outside the uncertainty margins.
Globally, extreme values of α should be avoided. For MUST, the best agreement with the reference partial similarity function is obtained with α 0.5, whereas α 0.66 seems to be better for IHOP. However the value of α 0.5 is a good compromise for the two cases, and makes it possible to have a resolved TKE closer to the theoretical curb in LES regime. Therefore, α 0.5 is retained as the default value in Eq. 7.

Cloudy Evolving Boundary Layer
In the previous section, α has been tuned on idealized cases of convective (IHOP, in Section 4.1) and neutral (MUST, in Section 4.2) boundary layers. In this section, L is tested on a cloudy evolving boundary layer.

Ability to Reproduce the Cloud Types and Large Scales
The CONSTRAIN case is more critical as it combines turbulence from thermal and dynamical origins, as well as presence of mixed clouds. CONSTRAIN was held in the North of Scotland. The simulation begins at 00 UTC -1 local hour. Figure 9 shows the cloud top and base temporal evolution in simulations undertaken with DEAR, RM17 and L at 125, 250, and 500 m grid spacing. Results are compared with the LES reference (in black) at 62.5 m grid spacing. The cloud base is well represented whatever the resolution and the mixing length. For all resolutions and simulation times, the boundary-layer height is well recovered although BL89 and RM17 result in a Frontiers in Earth Science | www.frontiersin.org January 2021 | Volume 8 | Article 582056 too shallow boundary layer at 125 m grid spacing (see also Figure 10). Figure 10 shows the vertical profiles of mean cloud fraction produced by the different mixing lengths in the stratocumulus stage (2 h local hour) and the cumulus stage (12 h local hour). One can see a good transition from stratocumulus to cumulus at 12 h for all simulations, with a decreasing cloud fraction and an increasing cloud depth (see Figure 9). Moreover, the different simulations produce a correct maximal cloud fraction in the stratocumulus stage despite a too low maximal cloud fraction at large scales (cf. Figure 11a). The maximal cloud fraction is also well-reproduced by L at fine resolution in the cumulus stage ( Figure 11b). However, for all the mixing length, almost half of the cloud fraction is missing at 500 m grid spacing at 12 h ( Figure 11f).
The main problem is that the coarser the resolution, the lower the cloud top ( Figure 9). This can be mainly explained by the fact that the simulations do not have any non-local turbulence scheme such as a mass-flux scheme, limiting the cloud development.
During the first hours of the LES simulation L MO is close to 300 m, the boundary-layer height is about 1.2 km and rolls appear in the simulations. The buoyancy flux is a quantity of particular interest in order to verify the schemes ability to reproduce the convective rolls of stratocumulus (at t 2 h, see Figure 10a) and cells of cumulus (at t 12 h, see Figure 10b). As explained in de Roode et al. (2019), the two minima of w ′ θ ′ v show negative values that will act as a sink term of TKE and impact the future development of both the stratocumulus and the cumulus layers. In the stratocumulus case, the cloud-base minimum is small and well-represented by L at all resolutions. The minimum of w ′ θ ′ v at the cloud top is gradually less negative as resolution decreases. In the cumulus case, the boundary-layer, the cloud layer and the cloud top are well-represented by L at all resolutions, especially the minimum at the cloud base showing the decoupling between the cloud and the surface. The largest variability of w ′ θ ′ v between the resolutions is obtained at the cloud top associated to the entrainment, but L at 500 m succeeds to reproduce this process. Figure 12 shows the resolved TKE produced by the different mixing lengths in the stratocumulus stage and the cumulus stage of CONSTRAIN. Focusing on the results of the coarse-graining LES (in black), it is observed that, similarly to the neutral case, the finer the resolution, the larger the resolved TKE. Contrary to the buoyancy where only a small sensitivity to the resolution has been FIGURE 8 | Sensitivity test to the α from 0.33 to 1 in IHOP (left) and MUST (right). The subgrid TKE is represented as a function of the resolution normalized by the ABL height. The reference partial similarity function for the TKE (Honnert et al., 2011) and the uncertainty margins are represented by the full and dashed black lines, respectively.

Resolved TKE Vs. Coarse-Grained LES
FIGURE 9 | Temporal evolution of the cloud base and cloud top in CONSTRAIN with DEAR, RM17 and L (red, green and blue lines, respectively) at 125, 250, and 500 m grid spacing (full, dashed and dotted lines, respectively).
Frontiers in Earth Science | www.frontiersin.org January 2021 | Volume 8 | Article 582056 observed, TKE is half less resolved at 500 m than at 125 m grid spacing. By construction, the L mixing length (in blue) depends on RM17 and DELT (close to DEAR). Figure 12 shows that the vertical profiles obtained with L are well located between those obtained with DEAR (in red) and RM17 (in green). More precisely, at high resolution (125 m grid spacing, Figures  12a,b) L mimics the behavior of DEAR when it mimics the behavior of BL89 and RM17 at 500 vm (Figures 12e,f).
In the stratocumulus stage (see Figures 12a,c,e), the L mixing length reproduces well the LES reference (in black) in the lower part of boundary layer. It is more difficult for this new mixing length to represent the stratocumulus cloud layer at 250 and 500 m, probably because neither L DEAR nor L RM17 are able to reproduce the characteristics of the stratocumulus at these resolutions. It should be noted that even in the cloudy part, L provides the profile closest to the reference. Be aware that we do not comment on the very first levels of the model, where the turbulence is mainly subgrid even with LES grid cells: as the turbulence is not mainly resolved, the LES cannot be considered as a good reference.
In the cumulus part (see Figures 12b,d,f), L reproduces well the resolved TKE profiles both in the cumulus layer and in the underlying boundary layer, particularly at fine scales. However, the cumulus present too large resolved turbulence at 500 m leading to an underestimated subgrid entrainment FIGURE 10 | Vertical profiles of mean cloud fraction for CONSTRAIN simulations at 2 h (left column) and 12 h (right column) at the grid spacing of 125 m (top row), 250 m (middle), 500 m (bottom row) with DEAR (in red), BL89 (in orange), RM17 (in green), L mixing length (in blue) and the reference LES coarse-graining (black).
Frontiers in Earth Science | www.frontiersin.org January 2021 | Volume 8 | Article 582056 ( Figure 10b). We assume that there is a lack of non-local turbulence for this resolution, whereas the boundary-layer thermals that produce the non-local turbulence are fully resolved at finer scales (Shin and Hong, 2014). This problem cannot be fixed by the mixing length itself (cf. Section 5).

DISCUSSION
As seen in this paper, L can be used in the LES regime. This is because, in practice, L and DEAR have the same behavior for all types of stability conditions. In unstable or neutral cases, DEAR is related to the grid size, which is equivalent to the use of the horizontal size of (LES nearly isotropic) grid meshes in L. In stable cases and in the inversion layers, DEAR limits the size of the eddies through a term accounting for the stratification. In L, the RM17 component plays this role (cf. Section 2.3). Over complex terrain, Cuxart (2015) underlines that hectometric scales using a 3D turbulence scheme are not in the LES regime, since the resolved motions fall outside the inertial sub-range and follow the definition of Smolarkiewicz and Prusa (2002) of the "Very-Large Eddy simulation" (VLES). In the same vein, Efstathiou et al. (2018) identified a regime called near gray zone, where most of the TKE is resolved but for which these coarse LES simulations should not be considered as LES converging because their grid length is not fine enough to present a clear inertial sub-range (see Sullivan and Patton, 2011). In a way, this paper deals with these VLES or near-gray zone regimes and L could be tested in cases of turbulence triggering at sun rise or turbulence in complex terrain.
The turbulence is not isotropic (as in LES) in the gray zone of turbulence, due to the residual traces of non-local turbulence. This is why Kitamura (2015) proposed two mixing lengths in the gray zone, one horizontal and one vertical. In the present paper, as in LES, the turbulence scheme is 3D isotropic at all resolutions by construction. L is built to be used in isotropic schemes, in order to be compatible with LES current schemes. In our case, the anisotropy of the turbulent fluxes is entirely driven by horizontal or vertical gradients of the atmospheric quantities.
At kilometric-scales, L behaves as L RM17 : it reproduces well the neutral boundary layer. However, in CBL simulations at resolutions coarser than 500 m, a non-local scheme remains necessary because L alone is not able to reproduce the CBL like the combination of turbulence and mass-flux schemes would. We believe that with the adequate non-local scheme, L can be used at all scales.
The value of the α tuning parameter has to be valid whatever the boundary layer stability as RM17 is supposed to adapt to the BL stability. However, when L is used in an other model than Meso-NH, α may be different than 0.5. Although, L is not very sensitive to the value of α around 0.5, the determination of the α value should be undertaken as any other tuning parameter of the given model.

CONCLUSION
This work seeks to improve the turbulence scheme at the gray zone of turbulence scales, particularly when the resolution is too coarse and the meshes too anisotropic for the LES configuration to give a good subgrid/resolved partitioning of the TKE. This article describes a new mixing length, L, built on a previously proposed approach to blend a LES and a mesoscale mixing length in order to produce the right amount of subgrid turbulence at hectometric scales: L min(αL Δ , L RM17 ), where α is a tuning parameter that may depend on the boundary layer characteristics (not on the resolution). The sensitivity test to α performed herein establishes α 0.5 as the default value. In this formulation, the mixing length is smaller than the smallest of the two components. The merged mixing length does not explicitly depend on the model resolution, but in practice it increases with increasing grid size up to the mesoscale, giving it a self-adaptation faculty. L is tested on three different boundary-layer cases and compared with LES and LES coarsegraining results. The first case, IHOP, represents a freeconvective boundary layer, with no significant impact from the wind. In this case, the gray zone of turbulence and global parameters such as the boundary-layer height or the resolved/ subgrid partition are well represented by the new mixing length. Secondly, MUST, a neutral case of purely dynamic turbulence is studied. In this case, the correct representation of subgrid/ resolved partition at all scales allows a good representation of the wind profile. Finally, in CONSTRAIN, a case of stratocumulus to cumulus transition in a cold-air outbreak, the turbulence is well represented inside the boundary layer by L in the gray zone of turbulence. This new isotropic turbulent mixing length can then be used for simulations of neutral and convective boundary layers ranging from the LES to meso-scale resolutions, including the VLES and gray-zone regimes.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation. Frontiers in Earth Science | www.frontiersin.org January 2021 | Volume 8 | Article 582056