A Non‐Dimensional Index for Characterizing the Transition of Turbulence Regimes in Stable Atmospheric Boundary Layers

The transition from moderate to weak turbulence regimes remains a grand challenge for stable boundary layer parameterizations in weather and climate models. In this study, a critical horizontal Froude number (≈0.28) is proposed to characterize such a transition, which corresponds to the development of quasi two‐dimensional pancake vortices. Traditionally defined stability parameters corresponding to the critical horizontal Froude number are estimated and are consistent with values in the literature. The critical horizontal Froude number can recover previously used height‐ and site‐dependent mean wind speed thresholds. These findings offer a way to constrain the validity range of Monin‐Obukhov similarity theory in numerical models for weather and pollutants dispersion.

• A critical horizontal Froude number (≈0.28) is proposed to characterize the transition from moderate to weak turbulence regimes • The critical horizontal Froude number corresponds to the development of pancake vortices • Previously used height-and site-dependent mean wind speed thresholds can be recovered from the critical horizontal Froude number

Supporting Information:
Supporting Information may be found in the online version of this article.
depressed and detached from the surface, leading to suppressed vertical mixing and the breakdown of bulk formulations and MOST.Knowledge on the transition between the two turbulence regimes is crucial for constraining the application range of MOST and for further improving turbulence parameterizations in stable boundary layers (Mahrt, 2014;Sandu et al., 2013).
Researchers have not agreed upon which index works best to characterize the transition between the two turbulence regimes in stable atmospheric boundary layers.The mean wind speed threshold (U s ), as a dimensional quantity, cannot be expected to be universal (Mahrt et al., 2015;Van de Wiel et al., 2012).The dual characteristics of heat flux and the stability parameter corresponding to the maximum heat flux (in terms of magnitude) (Lapworth & Osborne, 2020;Mahrt et al., 1998;Malhi, 1995), albeit widely studied and recognized, have been found to fail to separate turbulence regimes as the maximum heat flux tends to occur in the weakly stable regime (Acevedo et al., 2019).The Richardson numbers (including the gradient Richardson number and the flux Richardson number) have been extensively studied, motivated by the seminar work by Miles and Howard (Howard, 1961;Miles, 1961;Miles & Howard, 1964), yet the exact meaning of critical (or sometimes maximum) Richardson numbers remains to be clarified (Bou-Zeid et al., 2018;Galperin et al., 2007;Katul et al., 2014;Li et al., 2015).Other work proposed non-dimensional indices to characterize horizontal meandering (Mortarini et al., 2019) and vertical coupling (Peltola et al., 2021), but the exact transition thresholds for these indices remain empirically determined.This study aims to bridge this research gap by proposing a new non-dimensional index for characterizing the transition between moderate and weak turbulence regimes and determining the transition threshold value theoretically.In the following, the data set and methodology are elaborated in Section 2; a theoretical model is introduced in Section 3; interpretations of results are presented in Section 4; conclusions and discussion can be found in Section 5.

Data Set
The data from the Cooperative Atmosphere-Surface Exchange Study in October 1999 (CASES-99) field experiment are used (Poulos et al., 2002).The CASES-99 data set contains ultrasonic turbulence measurement on a 60-m tower in southeast Kansas (located at 37.649°N, 96.736°W) during the month of October 1999.The tower site was surrounded by the relatively flat terrain.Sonic anemometers were installed at 1. 5, 5, 10, 20, 30, 40, 50, and 55 m.Since the lowest sonic anemometer was moved on 20 October, the measurements at 1.5 m are not used in this study.Data are split into 15-min segments.

Data Processing
The coordinate system is defined that x, y, and z-coordinates represent the streamwise direction, the spanwise direction, and the vertical direction, respectively, and the corresponding velocity components are u, v, and w.Due to this rotation, the spanwise mean wind speed and the spanwise momentum flux become sufficiently small and are neglected.The mean wind speed is then defined as   =  , where the overbar represents the time average.Turbulent momentum fluxes and heat fluxes are computed as − ′ ′ and ′ ′ with the prime indicating the turbulent fluctuations from the time averages, and the friction velocity is calculated as , where θ denotes the potential temperature, g is gravity acceleration, and z denotes the distance to the ground.The mean wind speed gradient   =   and the mean potential temperature gradient    are obtained by fitting third-order polynomial functions to the mean wind speed and mean potential temperature profiles and then taking the derivative of the fitted functions.Using the above definitions, the Obukhov length is computed as The streamwise and vertical turbulence integral length scales are obtained by fitting exponential functions to the autocorrelation as (Salesky et al., 2013): where R uu (∆x,0) and R ww (0,∆z) denote the streamwise and vertical autocorrelation which are calculated as with α denoting u and w, ∆x, and ∆z denote the streamwise and vertical spatial lags, and  ℎ and   denote the streamwise and vertical integral length scales (Li, 2020;Salesky et al., 2013).An example of deriving   through fitting is given in Figure S1 in Supporting Information S1.It can be seen that the fitting method is reasonable for all reference heights except the lowest one (5 m), and thus the calculated   at 5 m is discarded.
The horizontal Froude number is defined as , where σ u is the standard deviation of streamwise velocity.
The buoyancy length scale is defined as , which represents the vertical displacement of a fluid parcel if all its kinetic energy were converted to potential energy (Billant & Chomaz, 2000b).In this sense, the horizontal Froude number can be regarded as the ratio of the buoyancy length scale to the streamwise turbulence integral length scale.Similarly, the vertical Froude number can be defined as , which represents the ratio of the buoyancy length scale to the vertical turbulence integral length scale.Note that the vertical Froude number is still defined with σ u , not the standard deviation of vertical velocity (σ w ).This is a major difference between our work and the study by Peltola et al. (2021).We choose σ u as again     represents the vertical displacement of a fluid parcel if all its kinetic energy were converted to potential energy (Billant & Chomaz, 2000b;Riley & Lindborg, 2008) while     only represents the vertical displacement of a fluid parcel within the buoyancy time scale.

Theory
Our starting point of characterizing the transition is to assume that the turbulence regime transition concurs with the development of quasi two-dimensional pancake vortices, which provides a constraint on the vertical length scale of turbulence.Vortex instability (Herring & Métais, 1989) in stably stratified fluid has been investigated by setting long vertical columnar vortex pair experimentally and numerically (Billant & Chomaz, 2000a, 2000b, 2000c;Leweke & Williamson, 1998).For weak stratification, the elliptic instability prevails by the gradual bending of each vortex core in the opposite direction to the vortex periphery.As stratification becomes stronger, the zigzag instability plays an increasingly important role.The vortex pair twists with almost no change of the dipole's cross-sectional structure and is ultimately sliced into thin horizontal layers of independent pancake dipoles.Under such conditions, turbulence behaves as quasi two-dimensional pancake vortices.However, we note that the dynamics of these pancake vortices remain closer to three-dimensional turbulence rather than exact two-dimensional turbulence (Sozza et al., 2015).
Since the shape of pancake vortices is buoyancy driven, the horizontal Froude number, as a stability parameter for assessing buoyancy effects, is often used to describe the development and evolution of pancake vortices with stratification (Mater & Venayagamoorthy, 2014a, 2014b).Further analysis showed that the typical thickness of pancake layers scales with the buoyancy scale (Billant & Chomaz, 2001;Lindborg, 2006;Riley & Lindborg, 2008).Thus, we expect F v approaches unity at the transition (i.e., F vc = 1 where the subscript c indicates critical) with the development of pancake vortices.However, it is much difficult to compute F v as it requires ultrasonic measurements at multiple heights and data fitting across heights (Figure S1 in Supporting Information S1).In contrast, the horizontal Froude number F h only requires ultrasonic measurements at a single height.Therefore, for practical purposes, it is better to formulate the transition in terms of F h instead of F v .
Building on this thinking, our model starts in the weak stratification and then identify the transition (i.e., a critical F h value) by introducing the constraint of F vc = 1 at the transition.By starting from weak stratification, many assumptions in our theoretical model should still be applicable.Essentially, we aim to derive a relation between F h and F v , which requires to derive a relation between L h and L v .This is accomplished by using a simplified spectral model introduced below.
While the Kolmogorov −5/3 power law scaling for isotropic turbulence ends at the Ozmidov scale that indicates the largest scale of isotropic turbulence (Li et al., 2016), Lindborg (2006) reported that a new −5/3 power law scaling holds for anisotropic turbulence bounded in pancake layers.The new −5/3 power law scaling in this so-called buoyancy subrange has received some support from atmospheric surface layer studies (Cheng et al., 2020), but the mechanism for the new −5/3 power law scaling differs from that in isotropic turbulence.Energy is first transferred directly from the large scales to the Ozmidov scale by the zigzag instability and Kelvin-Helmholtz instability, and then transferred to the Kolmogorov scale by the traditional energy cascade process described by the Kolmogorov theory (Deloncle et al., 2008).Given these results, we assume the −5/3 power law scaling extends from the integral length scale to the inertial subrange of isotropic turbulence.To facilitate analytical treatment, we further ignore the spectral drop in the dissipation range.This is justified given that the dissipation range contributes little to the total energy (variance).At the low frequency end, we assume that the spectra flatten once the scale becomes larger than the integral length scale.Although the −1 power law scaling is often observed in neutral conditions for the streamwise velocity spectrum (Calaf et al., 2013;Drobinski et al., 2004;Katul et al., 1995Katul et al., , 1996Katul et al., , 2012;;Marusic & Monty, 2019), previous studies also showed that such a −1 power law scaling ceases to exist under moderately stable stratification (Krug et al., 2019).To support this argument, the streamwise and vertical velocity spectra are computed from the CASES-99 data and are shown in Figure S2 in Supporting Information S1.Clearly, the low frequency components of both streamwise and vertical velocity spectra are closer to a zero power law scaling than a −1 power law scaling.Therefore, the streamwise and vertical velocity spectra are described as where E u (k) and E w (k) denote the streamwise and vertical velocity spectrum respectively,   = 18 55 × 1.5 and   = 24 55 × 1.5 are the Kolmogorov spectral constants for one-dimensional velocity spectrum, and k denotes the one-dimensional wavenumber.Integrating Equations 3 and 4 within 0 < k < ∞ yield Invoking the definitions for F h and F v gives To close Equation 7, we need to further examine the relation between  (Banerjee et al., 2016) when turbulence is separated from gravity waves.Namely, where σ v indicates the standard deviation of velocity in the spanwise direction.The 1:1 comparison between   2  and   2  +  2  is shown in Figure S3 in Supporting Information S1 to support Equation 8.The standard deviation of spanwise velocity in Equation 8 can be closed by the spanwise velocity variance budget for stationary and horizontally homogeneous surface layer flows which reads where ϵ v denotes the TKE dissipation rate in the spanwise direction, P denotes the pressure, and ρ denotes the air density.Vertical transport term is neglected because the magnitude of this term is small compared with other terms under mild to moderate stratification (e.g., Li et al., 2016).The term on the left-hand side of Equation 9 can be estimated as  1 3  because small-scale energy dissipation is isotropic (Lindborg, 2006).The term on the right-hand side of Equation 9 is the correlation between pressure fluctuations and velocity shears that describe how TKE is redistributed by pressure perturbations.Traditionally, Rotta's return-to-isotropy hypothesis (Rotta, 1951) is used to parameterize this pressure correlation term in atmospheric boundary layer modeling (e.g., Katul et al., 2014;Zilitinkevich et al., 2007), which is expressed as (Pope, 2000): where C r = 1.8 is the Rotta constant and τ is a relaxation time scale representing how fast TKE is dissipated.Here we formulate   = to focus on the spanwise TKE content relative to its dissipation rate.Equations 9 and 10 give the form of   2 as: Then, by substituting Equations 8 and 11 into Equation 7, the sought relation can be derived: At the transition,   = 1 , corresponding to a critical horizontal Froude number ℎ = 0.28. (13)

Results
Using the CASES-99 data, F v is shown as a function of F h in Figure 1.As one can see, Equation 12 is in good agreement with observations in weak stratification, supporting our theory.Note that a decrease in F h , which corresponds to moving toward the right of Figure 1, indicates an increase in the stable stratification.When F h < F hc , the calculated F v from Equation 12 deviates from the observed F v , which remains around unit.This is consistent with the key assumption invoked in our model, namely, the transition concurs with the development of pancake vortices in the stable atmospheric boundary layers.The mean deviation of the observed F v when F h > F hc is 0.35, implying that the uncertainties introduced by calculation of the potential temperature gradient and estimation of the L v have a relatively small impact on the results.To exclude the potential impact of self-correlation on Figure 1, the relations between L v , L b , and L h are shown in Figure S4 in Supporting Information S1.It can be seen that when F h > F hc , the relationship between L v and L h is reasonably described by Equation 12.When F h < F hc , L v is no longer correlated with L h and instead becomes equal to L b .In short, the observation evidence supports the use of a critical horizontal Froude number F hc = 0.28 to indicate the development of pancake vortices in stable atmospheric boundary layers.
Following Sun et al. (2016), the HOST patterns are shown in Figure 2a.The magnitude of u * is around 0.15 m s −1 or less in the weak turbulence regime and is above 0.15 m s −1 in the moderate turbulence regime.In this plot, the mean wind speed threshold used to separate the two turbulence regimes is height dependent.It is also site dependent as shown in other work (Mahrt et al., 2015;Van de Wiel et al., 2012).In contrast, when the u * is plotted against F h in Figure 2b, the critical horizontal Froude number F hc = 0.28 well separates the two turbulence regimes.More importantly, no height dependence is found.In this sense, F hc is a better, non-dimensional threshold separating the two turbulent regimes.
One might be tempted to characterize the transition using traditionally defined stability parameters such as z/L and Ri g .In the context of HOST patterns, Figures 2c and 2d shows that neither z/L nor Ri g is a good indicator given the large scatter.However, it is insightful to estimate the critical values of z/L and Ri g corresponding to F hc as z/L and Ri g are widely utilized in operational stable boundary layer parameterizations.To do so, all z/L and Ri g data points in the interval 0.26 < F h < 0.3 are considered, which give a mean value of 2.15 ± 0.39 for z/L (± refers to the 95% confidence interval) and a mean value of 0.241 ± 0.026 for Ri g .Note that we do not observe strong variations of these values using other F h intervals such as 0.27-0.29 or 0.25-0.31.These values, shown in Figures 2c and 2d, are consistent with those used in the literature (Galperin et al., 2007;Grachev et al., 2013;Li et al., 2015;Mahrt et al., 1998;Miles, 1961;Sorbjan, 2010;Stull, 1988;Sun et al., 2016).
Our results may shed some insights into the meaning of critical Richardson numbers often quoted in the literature.Miles (1961) reported that wave-like perturbations are dynamically stable when the flow is characterized by a linearized Euler equation and when Ri g exceeds some critical value.While this seemingly suggests that stable boundary layer flows would laminarize when Ri g exceeds some critical value, this association of critical Ri g with flow laminarization has been questioned (Galperin et al., 2007;Katul et al., 2014;Li et al., 2015;Zilitinkevich et al., 2007Zilitinkevich et al., , 2013)), namely, a critical Ri g associated with flow laminarization is not supported by observational and modeling evidence.Nonetheless, many studies did report a maximum (sometimes also called critical but is called maximum in our study to avoid confusion) flux Richardson number (Ri f ) that indicates substantial changes in flow behaviors.When the atmospheric stability is weak (e.g., when the Ri g is small), Ri f increases nearly linearly with Ri g .As the atmospheric stability becomes stronger, especially when the Ri g exceeds a value of about 0.2-0.25 (consistent with the values shown in Figure 2b), the Ri f no longer increases (or a maximum value for Ri f is attained).Previous studies have provided important theoretical insights into the meaning of this maximum Ri f .Using a reduced model based on the return to isotropy concept, Bou-Zeid et al. (2018) showed that the maximum Ri f is reached with zero vertical turbulent fluctuations (σ w = 0).Differently, Katul et al. (2014) and follow-up studies by Li et al. (2015) and Li (2019) interpreted the maximum Ri f as the threshold beyond which the Kolmogorov scaling is no longer a distinct feature for turbulent velocity and temperature spectra, which is consistent with observational data presented in Grachev et al. (2013).
In this regard, our results have three important implications.First, our results imply that the turbulent flow in stable boundary layers does not become laminar under strong stratification when Ri g exceeds 0.2-0.25 (see Figure 2d).Instead, the flow transitions to a weak turbulence regime characterized by pancake vortices, which provide a finite amount of turbulent mixing.Second, our model differs from that of Bou-Zeid et al. (2018) in that our model considers the structural characteristics of turbulence eddies and the transition identified by our model does not correspond to zero vertical turbulent fluctuations (i.e., two-dimensional turbulence).Previous studies have demonstrated that stably stratified turbulence still has a forward energy cascade, rather than an inverse energy cascade (Gucci et al., 2023;Sozza et al., 2015).Thus, it is important to recognize that although pancake vortices "look like" two dimensional, their dynamics remain three dimensional.Third, our model shares similarity with that of Katul et al. (2014) in that the spectra of horizontal and vertical velocity follow the Kolmogorov scaling at high wavenumber (Equations 3 and 4).The fact that our model results deviate from observational data when  ℎ < ℎ (Figure 1) suggests that  ℎ = ℎ is the point beyond which some key assumptions of our model break down.For example, velocity spectra in buoyancy subrange may differ from those described by Equations 3 and 4 when F h < F hc , as Kelvin-Helmholtz instabilities generated within pancake layers (Deloncle et al., 2008) can lead to the emergence of internal gravity waves.In this sense, F hc has similarities with the maximum Ri f suggested by Katul et al. (2014).However, we should emphasize that our model is formulated based on the horizontal Froude number while the model of Katul et al. (2014) yields an outcome regarding the flux Richardson number.Future work examining the connection between the critical horizontal Froude number, the maximum flux Richardson number, and internal gravity waves is welcome.
Lastly, we demonstrate that the mean wind speed threshold proposed by Sun et al. (2012) can be recovered from our theory.To do so we employ the stability-corrected logarithmic wind profile, which reads where the subscript s for u * and L indicates the transition at F hc = 0.28 (in practice the results at F hc = 0.28 are approximately estimated by the averages over the range 0.26 < F h < 0.3 to reduce uncertainty), a m = 5 is a constant, and z 0 is the roughness length which is taken to be 0.03 m following Van de Wiel et al. (2003).A comparison between the calculated U s and the mean wind speed threshold identified by Sun et al. (2012) (denoted observed U s ) is shown in Figure 3, with an average difference of the two profiles less than 0.5 m s −1 .In Equation 14, we assume u *s and L s to be height-independent, although in reality, u *s and L s are slightly height-dependent, which may contribute to the slight underestimation of U s below z = 30 m.

Conclusions and Discussion
The transition from moderate to weak turbulence regimes remains a challenge for stable boundary layer parameterizations in numerical weather and climate models.In this study, we propose a non-dimensional index for characterizing such transitions by invoking the development and evolution of pancake vortices.The length scale constraint of pancake vortices is introduced into a theoretical model, and the results show that a critical horizontal Froude number F hc = 0.28 is capable to characterize the transition of turbulence regimes, which is supported by the CASES-99 data set.The critical values of z/L and Ri g corresponding to F hc are also estimated in this study; the mean values with 95% confidence intervals are 2.15 ± 0.39 for z/L and 0.241 ± 0.026 for Ri g , respectively.
While a z-dependent mean wind speed threshold was utilized to identify the transition in previous studies (Sun et al., 2012), we find that the mean wind speed threshold can be recovered from F hc through the stability-corrected logarithmic wind profile.
The findings in this study have broad implications.For instance, Sandu et al. (2013) reported that the representation of stably stratified turbulence in numerical weather prediction models remains a challenge, especially for the weak turbulence regime.Ren et al. (2022) documented that the rapid accumulation of pollutants is related to weak turbulence during heavy haze.The knowledge of the transition of turbulence regimes can be directly used to constrain the validity range of Monin-Obukhov similarity theory in numerical simulations for weather and pollutants dispersion.Although numerical models do not directly provide values of L h , L h can be linked to σ u and ϵ through Equation 5, which are available in some turbulence parameterizations (Zilitinkevich et al., 2007(Zilitinkevich et al., , 2013)).
But future work needs to address how pancake vortices affect turbulent mixing efficiency and how to parameterize such mixing processes in numerical models.How pancake vortices differ from two-dimensional turbulence also needs to be studied.Moreover, the applicability of F hc over mountainous regions requires further investigations because large-scale internal gravity waves induced by topography can transfer energy to turbulence (Sun et al., 2015).This process can potentially disturb the slicing process by zigzag instability due to enhanced connection between adjacent layers.

Figure 1 .
Figure 1.The relation between vertical Froude number F v and horizontal Froude number F h .The red dots and vertical bars refer to the median values and corresponding standard deviations in F h bins, respectively.Each bin has a logarithmically equal width of 0.2.The blue solid line denotes Equation12.The horizontal and vertical black dashed lines indicate the unit F v and the critical horizontal Froude number F hc = 0.28, respectively.Note that x-axis is reversed so that moving toward the right of this figure indicates an increase in the stable stratification.

Figure 2 .
Figure 2. Relations between the friction velocity u * and (a) the mean wind speed U, (b) the horizontal Froude number F h , (c) z/L, and (d) the gradient Richardson number Ri g .Data points are medians in mean wind speed bins, which are the same as those shown in (a).The black dashed lines indicate critical values of (b) F h , (c) z/L, and (d) Ri g .Blue and red dashed lines in (c, d) indicate the upper and lower limits of the critical values at the 95% confidence, respectively.

Figure 3 .
Figure 3.The vertical profiles of calculated U s (blue line) by Equation 14.The red dots indicate observed U s in Sun et al. (2012).The shading indicates 95% confidence intervals for the calculated U s .