A Budget-Based Turbulence Length Scale Diagnostic

: The most frequently used boundary-layer turbulence parameterization in numerical weather prediction (NWP) models are turbulence kinetic energy (TKE) based-based schemes. However, these parameterizations suffer from a potential weakness, namely the strong dependence on an ad-hoc quantity, the so-called turbulence length scale. The physical interpretation of the turbulence length scale is difﬁcult and hence it cannot be directly related to measurements or large eddy simulation (LES) data. Consequently, formulations for the turbulence length scale in basically all TKE schemes are based on simpliﬁed assumptions and are model-dependent. A good reference for the independent evaluation of the turbulence length scale expression for NWP modeling is missing. Here we propose a new turbulence length scale diagnostic which can be used in the gray zone of turbulence without modifying the underlying TKE turbulence scheme. The new diagnostic is based on the TKE budget: The core idea is to encapsulate the sum of the molecular dissipation and the cross-scale TKE transfer into an effective dissipation, and associate it with the new turbulence length scale. This effective dissipation can then be calculated as a residuum in the TKE budget equation (for horizontal sub-domains of different sizes) using LES data. Estimation of the scale dependence of the diagnosed turbulence length scale using this novel method is presented for several idealized cases.


Introduction
The parametrization of turbulence processes in the Atmospheric Boundary Layer (ABL) is an essential component of numerical weather prediction (NWP) and global circulation (GC) models.
The parametrizations are based on broad physical understanding from observations, as well as on the statistical properties of resolved-turbulence flows from finer-scale model simulations. The grid spacing determines the extent to which the turbulent motions are resolved. In the traditional NWP and GC models all turbulent motions are sub-grid and turbulence needs to be fully parametrized. Large Eddy Simulation (LES) models have a smaller grid spacing and therefore, resolve a part of the turbulent motions, namely the large energy-containing anisotropic turbulent eddies. Therefore, only the smaller isotropic eddies need to be parametrized in LES. Between the range of grid-sizes of classical NWP/GC models and LES lies the gray zone of turbulence [1][2][3], where turbulence parametrization is more difficult, because the anisotropic turbulent eddies are only partly resolved and the other part of their impact needs to be parametrized. Until recently, there was no need for such a gray-zone representation of turbulence. However, increases in computational power and the assumption that higher resolution models produce better forecast push the NWP/GC models towards the need for a representation of the gray zone of turbulence. by [27], who use a prescribed length scale as a reference in a relaxation equation, which includes an advection and a turbulent transport term.
Separately from the above approaches, the turbulence length scale can also be interpreted from the perspective of the turbulence energy spectra (see e.g., [28]). The atmospheric turbulent kinetic energy spectra can be divided into three distinct sections: the energy production range, the inertial sub-range, and the dissipation range. In the energy production range, the TKE is generated by shear and generated or reduced by buoyancy and it is, therefore, dominated by anisotropic large eddies. In the isotropic inertial sub-range, the TKE is internally transported from larger eddies to smaller eddies. In the dissipation range, the TKE of the smallest eddies is dissipated by viscous effects. In the quasi-steady inertial sub-range, the net energy coming from the larger eddies is in equilibrium with the net energy cascading to smaller eddies, resulting in a constant slope of the log-log energy spectrum [29]. Reference [30] relate the turbulence length scale to a characteristic scale of the spectra. Some authors instead interpret the integral turbulence scale (characteristic size of the most energy-containing turbulence eddies) or the grid scale as the turbulence length scale [31,32].
Alternatively, the spectral form of TKE equation [4,33,34] can be used to estimate the turbulence length scale. The length scale is then usually associated with TKE loss via molecular dissipation and can be computed from the dissipation rate and the spectrally integrated TKE [15]. The dissipation rate can be determined from the constant slope of the log-log TKE spectra in the inertial sub-range [35]. The transfer of TKE across scales is though not considered, because it is either assumed that all turbulence motions are sub-grid, as in the case of low resolution NWP/GC regime, and hence the net sub-grid cross-scale transfer is equal to zero [28], or that the cross-scale TKE transfer is modeled explicitly, the case of the LES regime [36]. However, in the gray zone of turbulence, the effect of the sub-grid cross-scale TKE transfer should be accounted for and its influence could be included in the turbulence length scale.
The limitations of the TKE scheme in the gray zone can be partially solved by extending the TKE scheme with influence of the 3D turbulence [37,38] or by TKE blending [39]. Other authors propose to stay within the 1D TKE framework, and use a scale-dependent turbulence length scale. The turbulence scale is obtained by pragmatic blending between the NWP/GC length scale and a small-scale 3D mixing length scale [40][41][42]. The calibration of the blending is of critical importance in this approach.
In the present paper, we propose a new turbulence length scale diagnostic based on LES data that contains the influence of the sub-grid cross-scale TKE transfer. In the procedure, the sum of the molecular dissipation and the cross-scale TKE transfer are handled as an effective dissipation, which is associated with the new turbulence length scale. The effective dissipation is expressed as a residuum in the budget equations for the TKE and the scalar variances using LES data. The diagnosed turbulence length scale is thus resolution-aware and may serve as a basis for the development of turbulence length scale parameterizations for models operating in the gray zone of turbulence. This paper is organized as follows. Section 2 gives an overview of the spectral TKE equation and the formulation of the new turbulence length scale is derived. In Section 3 a LES-based diagnostic of the new turbulence length scale is proposed. The diagnostic is evaluated in Section 4 on LES data from several simulations. The results are discussed in Section 5.

Spectral TKE Equation
To understand the role of the molecular dissipation and the transfer of kinetic energy across scales, we look at the spectral representation of the TKE equation for horizontally homogeneous turbulence (see e.g., [4,34]): where k is the wave number, ν is the molecular viscosity, S x (k) are the cospectra of the corresponding covariances or TKE, τ r (k) is the convergence of TKE transport across the spectrum; and ADV and TT are spectral representations of the vertical advection and turbulent transport, respectively. The averaged quantities are marked by a bar and the fluctuations from the mean are marked by a prime symbol.
Integrating the TKE cospectra from some lower bound, k l , to infinity yields the sub-filter TKE: The cross-scale TKE transfer is obtained by integrating the convergence term (see e.g., upper panel of Figure 5.16 in [4]): The convergence term, τ r (k), is a gain term in the dissipation range, and is negligible in the inertial sub-range (LES regime). The viscosity term in Equation (1) is significant only in the dissipation range and can be neglected otherwise. In the energy production range, the buoyancy and shear terms are non-zero. In this regime, τ r (k) acts as a loss term and replaces the molecular dissipation as the main loss term in the spectral TKE equation [28].
For NWP/GC models all turbulence motions are unresolved. Thus, the sub-filter TKE is equal to the full TKE (including all turbulence scales): and the classical TKE equation for NWP/GC models is obtained by integrating the spectral TKE equation over the whole spectrum: where is the TKE dissipation rate: The cross-scale TKE transfer Tr(k l ) is equal to zero when computed over the whole turbulence spectrum: and equal to when integrated from a cutoff in the inertial sub-range, k i : For the gray zone of turbulence, the prognostic equation for the sub-filter TKE is obtained by integrating Equation (1) from a cutoff scale in the energy production range, represented by a wave number, k c : Compared to Equation (5), the molecular dissipation rate is replaced by an effective dissipation rate, c . The effective dissipation rate is then equal to the molecular dissipation rate minus the cross-scale TKE transfer T r (k c ) at the cutoff scale:˜ As described above, T r (k c ) is by definition equal to zero when all turbulence motions are sub-grid, or it is assumed to be modeled explicitly as part of the three-dimensional shear-production term in LES [36]. However, neither condition is met in the gray zone of turbulence for NWP/GC models. The cutoff wave number is in the energy-containing range, hence T r (k c ) is non-zero. Also, the one-dimensional (vertical direction) shear term parametrization in NWP/GC models is unable to fully model the complex three-dimensional interactions between the different spatial scales, responsible for the cross-scale TKE transfer [28]. Because T r (k c ) is always smaller than , the effective dissipation rate,˜ c , is always positive, representing a loss term in the TKE equation (−˜ c < 0).

Formulation of the New Turbulence Length Scale
Our goal is to replace the parametrization of the dissipation rate in a standard TKE scheme for NWP/GC models with a parametrization of the effective dissipation rate via modification of the turbulence length scale. The Reynolds averaged TKE equation for NWP (assuming horizontal homogeneity) that corresponds to Equation (5) is (see e.g., [4]): where x, y, and z are horizontal and vertical coordinates; u, v and w are the corresponding components of the wind velocity; t is the time, p is the atmospheric pressure, e k is the TKE: θ v is the virtual potential temperature; θ vr is a reference virtual potential temperature; ρ 0 is a reference density; g is the gravitational field strength.
The first term on the right-hand side in Equation (11) represents the vertical turbulent transport of TKE, including the pressure correlation term. The second term is the shear-production term, representing kinetic energy transfer from mean to the turbulent scales. The third term is the buoyancy production/destruction term.
Following [29], the dissipation rate in the prognostic TKE equation is parametrized with the help of the turbulence length scale or the corresponding time scale [9]: where L is the classical turbulence length scale, τ k is the dissipation time scale for TKE and C is a closure constant (see e.g., [10,43]).
In the gray zone, the size of the grid elements, l c , lies in the energy production range of the turbulence. Hence, the turbulent flow is not anymore entirely sub-grid, but is partly resolved and partly sub-grid. The resolved TKE is implicitly modeled via the dynamics and only the sub-grid (or sub-filter if block filter is used) scale TKE needs to be parametrized. To underline this fact, we denote the average over a grid element of size l c with the [] l c operator, and the fluctuations from this average with : where ψ stands for any relevant quantity.
In LES, the cross-scale TKE transfer from the resolved scales to sub-grid scales is modeled together with the kinetic energy transfer from mean to the sub-grid scales via the three-dimensional gradient production term [36]. However, we intend to use the NWP/GC framework, where parametrizations are computed in vertical columns without any unresolved horizontal transport. Thus, the shear-production term for the sub-grid scale TKE is one-dimensional and contains only part of the cross-scale TKE transfer. To account for the influence of the cross-scale transfer in the sub-grid scale TKE equation, the molecular dissipation rate, , needs to be replaced with a reduced (by the one-dimensional resolved part of the cross-scale transfer) effective dissipation rate,˜ c (l c ) (see Equation (10)).
For this purpose, we define a new turbulence length scale, valid also in the gray zone. In analogy with the traditional turbulence length scale, L , defined by: the new turbulence length scale, L c , depends on the sub-gridscale TKE, e l c k , and the effective dissipation rate,˜ c (l c ): where e l c k is the sub-grid scale TKE corresponding to a grid of size l c . The corresponding sub-grid scale TKE equation becomes: The advantage of Equations (16) and (18) is that they account for the cross-scale TKE transfer (as is the case in Equation (9), but formally do not modify the classical TKE equation for NWP/GC models (Equation (11)).

Estimation of the New Turbulence Length Scale
We estimate the turbulence length scale influenced by cross-scale TKE transfer, L c , using LES data. The traditional turbulence length scale, L , can be diagnosed for different scales from LES data for a specific NWP/GC model resolutions according to Equation (15): where: Here e k (k c ) is the TKE obtained by integrating the TKE cospectra from the smallest scales to the NWP cutoff. e k,i is the TKE corresponding to wave numbers larger than k i (cutoff wave number in the inertial sub-range, see Equation (8), and is equal to the sub-grid scale TKE in the LES model. The dissipation rate, , is available directly from the LES data and the TKE cospectra can be computed from the velocity fields via Fourier transformation.
The new scale-dependent length scale, L c , (see Equation (16)) can in principle be estimated in a similar way: It is, however, costly to compute the contribution of the cross-scale TKE transfer to the effective dissipation rate in spectral space. Instead of integrating the TKE equation in spectral space (application of a clean wave-cutoff filter), we use volume averaging (application of a box filter in physical space) to compute the terms in the prognostic TKE equation (Equation (18)) and express the reduced effective dissipation rate,˜ l c , as the residuum (see e.g., [44]): The terms on the right-hand side of Equation (22) are computed from the LES data. To focus on the dependence of the length scale on horizontal resolution, the [] l c operator is interpreted as averaging over a rectangular horizontal sub-domain of size l c on each vertical model level.
Thus, the new scale-dependent length scale, L c , is computed by averaging over all sub-domains with size equal to l c (similar to [15]; see Figure 1 for TKE example): where <> denotes the averaging operator over all sub-domains with size l c . For comparison, the classical turbulence length scale is estimated in the same manner: where is assumed to be equal to the effective dissipation rate,˜ c , estimated from Equation (22) for the whole domain. The number and the size of the sub-domains are chosen so that both samples are large enough to ensure the representativeness of the resulting quantities.
Also, LES domain size is chosen large enough to include a sufficient number of sub-domains with spatial scales above the energy production range of TKE. The smallest sub-domains have a spatial scale in the gray zone of turbulence. In this way, the scale dependence of the diagnosed turbulence length scale can be investigated from conventional to gray-zone resolutions.
To verify whether the sub-domains reach into the energy production range, their sizes are compared with the longitudinal integral turbulence length scales (see e.g., [45]): which is an estimate of the size of the most energetic turbulent eddies. The integral length scales are computed from data of the whole domain, and the integration is done only up to the value where the auto-correlation function falls off to 1/e [45].

LES Simulations
The new length scale diagnostic is tested on four idealized cases, all developed by the Global Energy and Water Cycle Experiment (GEWEX): the drizzling stratocumulus case based on the first research flight (RF01) of the second Dynamics and Chemistry of Marine Stratocumulus (DYCOMS-II) field study [46], the continental cumulus case based on measurements over the Atmospheric Radiation Measurement (ARM) program Cloud and Radiation Testbed (CART) site in Oklahoma [47,48], the trade wind cumulus case based on observations from the Barbados Oceanographic and Meteorological Experiment (BOMEX) [49], and the fourth test case is based on the GEWEX Atmospheric Boundary Layer Study (GABLS1) project [50][51][52].
This relatively diverse set of cases is chosen to demonstrate the capabilities of the new turbulence length scale diagnostic in different ABL regimes. The results can additionally indicate, to what extent the turbulence length scale is case sensitive.
All simulations are performed with the MicroHH LES model [53,54]. MicroHH is setup and run according to the description of the cases. Moist thermodynamics is employed in the simulations for DYCOMS-II, ARM and BOMEX case. The GABLS1 simulation is run with dry thermodynamics. Parametrization of radiation and microphysics is used only in the DYCOMS-II case. These processes are deactivated in the remaining cases.
With the aim of gaining statistics for a wide range of horizontal scales, a relatively large horizontal domain size and high horizontal resolution are used for the LES runs. This enables study of the influence of the horizontal scale on the diagnosed turbulence length scale. The vertical domain size and resolution follow the description and recommendation in their respective inter-comparison studies (while maintaining a physically reasonable aspect ratio of the gird box). The resolutions, the sizes of the domain and the integration time are summarized in Table 1. The turbulence length scale is diagnosed on data from rectangular horizontal sub-domains. We construct the sub-domains by decomposing the whole domain into equal squares, with 1, 2, 4, 8, 16, and 32 sub-domains in each direction.
Example of sub-domains construction is displayed in Figure 2.
For the turbulence length scale diagnostic, the TKE budget terms are averaged over all the sub-domains of the same size. The LES statistics are averaged over one hour in time.
All the MicroHH LES simulations reproduce the results of the respective LES inter-comparison studies.

Fit of Turbulence Length Scale
For demonstration, the diagnosed turbulence length scales are fitted with the following algebraic height-dependent relationship [16,18] : where κ = 0.4 is the von Kármán constant, C K = 0.0882 and C = 0.871 are closure constants, H abl is the ABL height; and a m , b m , β m and λ m are shape constants. The length scales are fitted with the least-squares method implemented in the scipy Python package. H abl and the shape constants are used as fitting constants. The bounds and the value of the first guess for the fitting constants are presented in Table 2. Table 2. First guess and bounds for the fitting constants in Equation (26).

Computation of Turbulence Fluxes
For illustration, the diagnosed turbulence length scales are used to compute the turbulence fluxes for the GABLS1 case, following the local down-gradient NWP turbulence parametrization of [11]: where C 3 is a closure constant. The flux Richardson number, Ri f , the TKE, the potential temperature θ, and the wind components are obtained by time averaging the corresponding LES data. The turbulence length scale L stands for either the diagnosed turbulence length scale, L c , or the parametrized turbulence length scale, L B , according to Equation (18) in [16]. L B is computed as a minimum of an algebraic [18] and a TKE-based [25] turbulence length scale relationship: where L BL is given by the balance between TKE and local stability. L up and L down are upward and downward free paths, respectively.

Gravity Waves
The velocity fluctuations in the LES data are caused by both turbulent motions and by gravity waves [15]. The effect of gravity waves is particularly strong above the ABL in the convective cases (BOMEX and ARM). Since it is difficult to separate the TKE from the wave kinetic energy in this region, the estimation of TKE from the LES data is inaccurate. Therefore, we avoid this problematic region in the fitting procedure, and perform the fitting only on data below the cloud top. Cloud top and cloud base are defined in this paper as layers above/under the cloud, where the cloud fraction becomes smaller than 0.005.

Results
Results of the turbulence length scale diagnostic for the four idealized cases are presented next. We start with the input parameters to the diagnostic: the TKE and the effective dissipation rate (see Equation (23)). The vertical profiles of the sub-gridscale TKE, e l c k , and the normalized differences in the effective dissipation rate between the sub-domain and the whole domain, ∆˜ c (l c ) / max = (˜ c (l c ) − )/ max z ( ), are presented in Figure 1. The first column in Figure 1 also shows the TKE budget terms for the smallest sub-domain. The changes in the size of the sub-domain are interpreted as changes in the corresponding horizontal grid resolution in an NWP model (see Section 3).
The TKE dissipation rate, , is set to the value of the effective dissipation rate of the whole domain, whose scale is assumed to be above the energy-containing range. The LES data is sampled every 60 s. The LES statistics are obtained by averaging over the last hour of the integration for the ARM, BOMEX, and DYCOMS-II cases. For GABLS1, the LES statistics are obtained by averaging over 45 min after 5 h of integration. To avoid division by zero in the following length scale diagnostic,˜ c is limited by a minimal value of 5 × 10 −5 m 2 s −3 , much smaller than the typical values of˜ c .
We can see in Figure 1 that the TKE, e l c k , gradually decreases with the decreasing size of the sub-domain throughout the entire vertical profile. In comparison, there is more vertical variability in the changes in the effective dissipation rate, because these are dictated by the budget of the TKE. Higher values of source terms and transport terms are usually correlated with stronger changes in˜ c . Among the four cases, the GABLS1 case manifests the weakest scale dependence due to relatively low intensity of the TKE source terms. The remaining cases show a similar vertical pattern with increasing horizontal grid resolution for an NWP model: increase of˜ c near the surface layer, strong decrease of˜ c under the cloud base and under the cloud top, and low scale dependence around the middle of the cloud layer.
The longitudinal integral length scales (see Equation (25)), L u,x i and L v,y i , the diagnosed turbulence length scale, L c , its fit according to Equation (26), L f it , and the parametrized turbulence length scale according to Equation (31), L B , are presented in Figure 3. The scale dependence of the classical length scale, diagnosed from , L (l c ) (see Equation (24), black and gray lines in Figure 3), is determined solely by changes in TKE. Therefore, the changes in L mirror the changes in TKE and are thus relatively uniform with height. L is, as expected, a monotonous function of the sub-domain size that converges linearly with e 3/2 k to zero (see Equation (24)) for very small sub-domains. The behavior of the new length scale L c (l c ) (see Equation (23)) is more complicated, because it depends on two factors, the TKE and˜ c . With an increase in resolution, the decrease of˜ c increases L c in comparison with L . This effect can even lead to a growth of L c with increasing resolution. The dependence of L c on the sub-domain size is thus not always monotonous and can result in local maxima of L c for a specific resolution within the gray zone.
Of course, L c should ultimately converges to zero for very small sub-domains. Reflecting the changes of˜ c , the resolution-dependence of L c changes with height. Obviously, the biggest differences between L c and L are in the regions where the magnitude of ∆˜ c is large: just below the cloud base and the cloud top in the three cloudy cases.
The changes in˜ c start to have a significant impact on the turbulence length scale diagnostic from a sub-domain size of 1.6 km × 1.6 km for the cloudy cases and at a sub-domain size of 50 m × 50 m for GABLS1. This can be interpreted as the beginning of the gray zone of turbulence. The largest impact of c changes are to be expected at sub-domain sizes around the longitudinal integral length scale, which is approximately 200 m to 400 m for ARM and BOMEX, 100 m to 300 m for DYCOMS-II and 10 m to 30 m for GABLS1. For all cases, the smallest sub-domain sizes are slightly larger than the integral length scales. Hence, all sub-domain sizes are in the upper part of the energy-containing range of TKE, or above it. The setup of the sub-domain is thus suitable for the scale dependence analysis of the upper part of the gray zone.
The influence of˜ c makes L c more variable than L , resulting sometimes in atypical vertical profiles compared to a classical turbulence length scale formulation. This is demonstrated via the fitting of L c with an algebraic height-dependent relationship (see Equation (26)). The fit works very well for the GABLS1 case, where L c and L are very similar. In the cloudy cases, the curvature of the vertical profile of L c increases with resolution, especially under the cloud base. This strong curvature cannot be completely fitted with the prescribed algebraic formulation.
Additionally, it would be impossible to fit a second peak in the vertical profile of L c with the chosen algebraic formulation, because the formulation implicitly implies a single-peaked vertical distribution. However, the existence of a second peek under the cloud top in our data is disputable. The presence of gravity waves near the top of the ABL in ARM and BOMEX, make the estimation of TKE and˜ c unreliable there (see Section 3.5). The proper resolution of this issue would be to separate the velocity fluctuations caused by gravity waves from the turbulence fluctuations and diagnose L c arising only from the turbulence part of the flow. This is however beyond the scope of this paper. Therefore, we opt for a pragmatic approach, and we fit only the data under the cloud top for the ARM and BOMEX case. Concerning the DYCOMS-II case, it is well known that most LES models have difficulties in modeling stratocumulus cases, especially near cloud top [55,56]. This is also seen in our simulations. Hence the TKE budget analyses for the DYCOMS-II case near the cloud top is also not reliable. Data from runs with an enhanced LES code could improve the diagnostic of L c for this case.
The estimation of both L c and L near the cloud top is strongly affected by inaccuracies in the data, because in this region the turbulence length scale is estimated as the ratio of two small numbers. To avoid this problem, the effective dissipation rate or the dissipation rate could be fitted/parametrized in a different way (e.g., direct parametrization without the need for a length or time scale ). This would, however, require a modification of the classical TKE scheme, which is beyond the scope of this paper.
To verify whether the diagnosed turbulence length scale is usable in a NWP turbulence scheme, it is compared with the parametrized turbulence length scale, L B (see Equation (31)). The diagnosed turbulence length scale for the larger sub-domains has a similar magnitude as L B . This confirms the validity of the turbulence length scale diagnostic for the NWP scales. There are differences in the shape between L c and L B , especially for the GABLS1 case. It is expected that the performance of an NWP turbulence parametrization would improve if L c or its fit would be used instead of L B for the computation of the turbulence fluxes (see Equations (27)-(29)). It must be however noted that the turbulence length scale is only one factor in the computation of the turbulence fluxes. The calibration of the closure constants, the shape of the stability functions, the choice of the stability parameter, and the parametrization of non-local effects also influence the quality of the turbulence scheme. Hence, the turbulence scheme should be tested as a whole and should be adjusted if a new turbulence length scale formulation is used. Such modifications are, however, beyond the scope of this paper.
For illustration, and instead of testing the fit of L c in an NWP parametrization, L c is used for the computation of the local down-gradient turbulence fluxes according to Equations (27)-(29) directly from the LES data. These fluxes are then compared with the turbulence fluxes computed using L B and with fluxes obtained directly from LES (see Figure 4) for the GABLS1 case. The cloudy cases are no included, because they would involve an additional parametrization of the non-local turbulence fluxes. The comparison in Figure 4 shows that the fit of L c and L outperforms the parametrized L B in GABLS1. The turbulence fluxes computed with L c are in good agreement with the LES fluxes and show a similar scale dependence. As expected, the computed turbulence fluxes are getting smaller for smaller sub-domains, because both the TKE and the turbulence length scale are getting smaller. From the perspective of the NWP parameterization, this corresponds to the increase of the horizontal resolution, where larger turbulence eddies are being resolved. To avoid double counting, the computed sub-grid scale turbulent fluxes should become smaller. The difference between L c and L for the smaller sub-domain (for larger sub-domains they are basically identical) shows that for L the decrease of the turbulence fluxes is stronger than for L c . This is the direct consequence of accounting for the cross-scale TKE transfer in the turbulence length scale diagnostic, which results in a smaller effective dissipation rate weaker in comparison to the dissipation rate.
When compared with the LES results, the scale dependence of the turbulent fluxes is too strong for both L and L c . However, the scale dependence caused by L c is smaller than the scale dependence cause by L and is thus closer to the LES results. Quantification of scale dependence is important for future development of a scale-dependent turbulent length scale formulation.  Figure 1, but comparison of the longitudinal integral turbulence length scales (see Equation (25)) for the zonal wind, L u,x i , and the meridional wind, L v,y i is computed for the period one hour before the end of the simulation (a,d,g,j); diagnosed turbulence length scales computed from the effective dissipation rate, L c (l c ) (see Equation (23)); and the dissipation rate, L (l c ) (see Equation (24)), (b,e,h,k); and their respective fits with the relationship in Equation (26) (c,f,i,l). The parametrized turbulence length scale according to Equation (18) in [16], L B , is also presented (c,f,i,l). The fits for the ARM and BOMEX cases are done only with data below the cloud top (indicated by horizontal line).

Discussion
We have derived and presented a new turbulence length scale diagnostic based on LES data that contains the influence of the sub-grid (in an NWP context) cross-scale TKE transfer. Differently from the classical turbulence length scale, L [15], (see Equation (15)) the new turbulence length scale, L c , is computed from the effective TKE dissipation rate,˜ c (l c ), instead of the TKE dissipation rate, . The effective dissipation rate takes into account the influence of the cross-scale transfer of TKE from the resolved to the sub-grid scales, which is non-zero in the gray zone of turbulence.
The effective dissipation rate is estimated from the TKE budget (see Section 2). The diagnostic of the associated turbulence length scale, L c , thus aims directly at the role of the turbulence length scale: parametrization of the TKE loss via dissipation and cross-scale TKE transfer. Thanks to this, L c does not require other simplified, limiting assumptions, compared to other formulations of the turbulence length scale (see Section 1). Consequently, our new approach offers a better reference for the parametrization of the turbulence length scale in the gray zone of turbulence, where neither the classical NWP or the LES assumptions are applicable. The new diagnostic can thus be used as a basis for a new scale-dependent turbulence length scale formulation suitable for the gray zone of turbulence.
The influence of the cross-scale TKE transfer in the upper part of the gray zone of turbulence on˜ c is presented in Section 4 for four LES cases: ARM, BOMEX, DYCOMS-II, and GABLS1. The differences between˜ c and basically indicate the largest scales at which the gray zone of turbulence starts to be significant, which in itself could be useful for gray zone detection.
Results in Section 4 show the shape of the vertical profile of L c varies more strongly with the horizontal scale than the shape of the L profile, due to the influence of˜ c on L c .
In the gray zone, the vertical profiles of L c have thus a rather atypical shape and cannot be fitted with the usual algebraic height-dependent relationships.
To fit the diagnosed L c values a modification of the algebraic formulations would be required, such as an increase in curvature and possibly the addition of a secondary maxima. We plan to propose a new height-dependent algebraic formulation based on our results in future research. We must however point out that finding a generally valid algebraic formulation, might prove to be difficult, because the properties of L c are case sensitive. Therefore, this will require a considerable effort and the analyses of more simulations. Another possibility is to use a TKE dependent, or prognostic turbulence length scale formulations [21][22][23][24][25]27] (see Section 1), which correspond better to the actual state of the turbulent flow.
In the gray zone, the turbulence length scale should be a function of the grid resolution. One approach is the blending of an NWP/GC turbulence length scale and a small-scale 3D mixing length scale (see Section 1) [40][41][42]. These methods deliver a turbulence length scale that is a monotonous function of the resolution. The blending is also a height independent procedure. We have shown in Section 4 that the classical turbulence length scale, L , is a monotonous function of the resolution and that its scale dependence does not change with height. These two properties are not always valid in the gray zone, as demonstrated by our new turbulence length scale, L c . This finding should be considered in the design of a new turbulence length scale parametrization in the gray zone of turbulence.
Our method of diagnosing the turbulence length scale for the purpose of NWP parameterization is universal, and should work for any turbulent flow. Nonetheless, the quality of the method is affected by the accuracy of the LES data. Above all, sufficiently large samples need to be used in the computation of the TKE budget terms. For this reason, relatively large horizontal domains and high horizontal resolutions are needed for the LES runs. Reliability of the diagnostic can be additionally affected by gravity waves, as has been discussed in Section 4. Separation of turbulence and gravity waves at the same scale is difficult, but, if possible, this would address the problem of the diagnostic of the turbulence length scale near cloud top for shallow convective cases. It should also be mentioned that LES models have their own inherent limitations. For example modeling of stratocumulus cases is not satisfactory with most LES models. There are proposals to mitigate this problem [55,56], and we would like to test data produced with such an enhanced LES code.
We have demonstrated that the fit of the diagnosed turbulence length scale has similar magnitudes as the parametrized turbulence length scale according to [16] (see Figure 3) and it performs well in a local down-gradient turbulence flux computation (see Figure 4). However, to assess the added value of the new turbulence length scale formulations, it needs to be tested in a full NWP parametrization. Acknowledgments: The authors would like to thank Chiel C. van Heerwaarden and T. J. Anurose who helped as technical support with running the LES code. The data for this study were generated with the LES model MicroHH, openly available in Zenodo at https://doi.org/10.5281/zenodo.822842 [54]. The configuration files and output statistics of the four simulations (ARM, BOMEX, GABLS1, and DYCMS-II) are openly available in Zenodo at https: //doi.org/10.5281/zenodo.3685801 [57]. The simulations were performed on MISTRAL, DKRZ, Hamburg, Germany.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript: