Detecting Nonequilibrium States in Atmospheric Turbulence

: In this work we show how to retrieve information about temporal changes of turbulence in the atmosphere based on in situ wind velocity measurements performed by aircraft. We focus on the stratocumulus-topped boundary layer high-resolution data taken by a helicopter-borne platform Airborne Cloud Turbulence Observation System (ACTOS). We calculate two nondimensional indicators, the dissipation factor and the integral-to-Taylor-scale ratio, and study their dependence on the Taylor-scale-based Reynolds number. By analyzing these results, we can identify regions where turbulence is in its stationary state and regions where turbulence decays in time or, on the contrary, becomes stronger. We can also detect nonequilibrium turbulence states, which indicate the presence of rapidly changing external conditions. SIGNIFICANCE STATEMENT: The purpose of this work is to retrieve new information on turbulence in the atmosphere. We consider data from a ﬁ eld study based on in situ observations from a measurement payload below a slow-ﬂ ying helicopter. We show that it is possible to retrieve information on temporal tendencies in a given region from such measurements. We attempt to estimate whether turbulence is in its stationary state, or possibly, turbulence decays due to insuf ﬁ cient production, or, on the contrary, becomes stronger. Such information provide useful knowledge of small-scale processes in the atmospheric boundary layer and can be used to improve their parameterizations.


Introduction
Atmospheric turbulence is a complex phenomenon, characterized by the presence of a plethora of scales (eddies) ranging from small viscosity-dominated to large, energy-containing structures.Turbulence may undergo space and time variations due to rapidly changing external conditions, it may be locally suppressed or enhanced.To describe characteristic features of turbulence, statistical theories are sought for.In this context, a number of recent research works (Vassilicos 2015;Mazellier and Vassilicos 2008;Bos and Rubinstein 2017;Bos et al. 2007;McComb et al. 2010) address the problem of the equilibrium Taylor law (Taylor 1935) and its failure in the presence of rapid changes of the system.A new, nonclassical, although universal scaling is introduced to describe the latter.
The Taylor law is a cornerstone relation which connects the statistics influenced by small scales, that is, turbulence kinetic energy dissipation rate with the statistics determined mostly by large turbulence eddies: the root-mean-square of the turbulent velocity fluctuations U and the integral length scale L where 5 2n s ij s ij and s ij 5 (u i /x j 1 u j /x i )/2, U 2 5 1/3 u i u i 5 2/3k, u i is the velocity fluctuation component, k is the turbulence kinetic energy, and n is the kinematic viscosity.
The assumption C 5 const underlies the well-known Kolmogorov-Richardson equilibrium cascade picture.It is also a foundation of many eddy-viscosity turbulence models, such as k- (Jones and Launder 1972), where the characteristic length scale of turbulence eddies is estimated from the turbulence kinetic energy k and the dissipation rate .Another conclusion that follows from Eq. ( 1) is that where k 5 15n/ √ U is the Taylor length scale and Re k 5 Uk/n is the Taylor-scale-based Reynolds number.Equation (2) implies that the spectral range between the small and the large eddies increases with Reynolds number, which is a wellknown property of turbulence.
The general validity of the "equilibrium" equations, Eqs. ( 1) and ( 2), was coming into question, as laboratory studies by Antonia and Pearson (2000) and Burattini et al. (2005) revealed that C varies considerably.This behavior was first explained by the dependence on the inflow and boun dary conditions.In another work by Bos et al. (2007) it was shown based on theoretical arguments that C takes different values in forced and decaying turbulence, even though the classical, Kolmogorov 25/3 form of the spectrum was assumed.
Next, breakthrough research works (Seoud and Vassilicos 2007;Valente andVassilicos 2011, 2012) revealed an "unusual," although universal dissipation scaling in laboratory experiments of decaying grid turbulence.Those authors argued that C is not constant, but depends on the inlet conditions and local Reynolds number where Re 0 5 U ' L b /n, U ' is the flow speed at the inlet, L b is the length scale defined by the grid spacing, m and n are power coefficients and it was found experimentally that both are close to unity.Moreover, instead of (2), another, striking relation was observed: that is, the spectral range between the large-scale and the viscosity-dominated region remained constant, even though Re k increased or decreased.Rubinstein and Clark (2017) argued that turbulence characterized by the classical laws ( 1) and ( 2) should be called "equilibrium" turbulence.In such type of flow the energy spectrum takes a self-similar form and the classical 25/3 scaling holds (Kolmogorov 1941).Nonequilibrium states of a flow field appear after a sudden change of external conditions (e.g., forcing) when the system evolves toward another equilibrium (Mahrt and Bou-Zeid 2020;Wacławczyk 2021Wacławczyk , 2022)).In the nonequilibrium turbulence, the turbulence dissipation coefficient C and the integral-to-Taylor-scale ratios are described by relations (3) and (4) and a nonequilibrium correction should be added to the spectrum.Bos and Rubinstein (2017) derived the form of this correction, as well as the nonequilibrium scaling laws for C and L/k based on theoretical arguments.
The nonequilibrium relations (3) and ( 4) are substantially different from their equilibrium counterparts in (1) and (2).Hence, classification of turbulence in a given area can be made on this basis.Previous research works which focus on such analysis concern laboratory or numerical experiments (e.g., decaying grid turbulence, turbulent boundary layer) under controlled conditions (Valente and Vassilicos 2015;Obligado et al. 2016;Cafiero and Vassilicos 2019;Obligado et al. 2022).The main idea of this work is to use this procedure to analyze data from in situ airborne measurements of atmospheric turbulence, where the flow conditions in the investigated region are unknown a priori and only limited information, that is 1D intersections of the flow field along the flight track are available.In principle, such analysis could be performed on any high-resolution observations of the atmospheric boundary layer.In this study we emphasize theoretical aspects of the method to be used; hence, we decided to analyze data already known from previous campaigns, but in a way that adds one more important aspect of turbulence.Our choice is the stratocumulus-topped boundary layer (STBL).Several research campaigns and subsequent data analyses were devoted to the study of the turbulence in STBL (see, e.g., Stevens et al. 2003;Malinowski et al. 2013;Siebert et al. 2010;Jen-La Plante et al. 2016).We focus on the data from the Azores Stratocumulus Measurements of Radiation, Turbulence and Aerosols (ACORES) campaign collected in the eastern North Atlantic (Siebert et al. 2021).Nowak et al. (2021) compared the properties of turbulence between the cases of coupled and decoupled STBLs [see Wood et al. (2015) for the definition of decoupling].They suggested that the observed differences in inertial range scaling might be related to different stages of turbulence lifetime (e.g., development and decay) but did not study these mechanisms in detail.In this work, from the same dataset, we derive all turbulence statistics necessary to calculate C , L/k, and Re k .The aim is to divide the flow area into regions where turbulence decays, develops or is in a stationary state, by investigating values of C .Moreover, by plotting C and L/k as a function of Re k we want to identify whether turbulence is in or out of equilibrium.Additionally, by the curve fitting of the data, estimation of the initial (or equilibrium) Reynolds number becomes possible.We also propose a procedure to roughly estimate the characteristic time scale of the turbulence decay/development.To the best of the authors' knowledge, this is the first application of such analysis in the study of atmospheric turbulence.
We perform a sensitivity study to show that all statistics, including , can be estimated from the wind velocity records, even in the case of nonequilibrium, provided that the resolution of the data is good enough and the flight segments are sufficiently long.These conditions are fulfilled by the ACORES data.Moreover, we extend analysis of Bos (2020) and estimate the upper bound of C for flows with nonzero buoyancy.
The present paper is structured as follows.In the following section 2 formal definitions of the "equilibrium" and "nonequilibrium" turbulence are presented and the theoretical derivations of Bos and Rubinstein (2017) are recalled.In section 3 the investigated data are briefly presented.This is followed by the detailed description of the methods applied in the data analysis.Results are presented in section 4, and finally, conclusions and perspectives are discussed in section 5.

Theory of nonequilibrium turbulence a. Nonequilibrium scaling laws
This section is devoted to the theory presented previously in the works of Bos and Rubinstein (2017) and Rubinstein and Clark (2017).Their derivations are briefly recalled here for the sake of clarity.
According to Rubinstein and Clark (2017) the term "equilibrium" turbulence, with C 5 const refers to the case where turbulence energy spectrum takes a (possibly timedependent) self-similar form: where k is the wavenumber, L(t) is a characteristic length scale (e.g., integral), and F is a function of a "fixed" form.All considerations are performed under the assumption of isotropy of the flow field.With the use of a spectral model those authors conclude the Kolmogorov K41 25/3 scaling is the solution for a stationary flow.In the case of nonstationarity, a small, subdominant correction appears; however, the K41 inertial range can still be treated as an approximate solution.Equation (5) hence describes turbulence in a state of a possibly unsteady equilibrium.As shown by Bos et al. (2007) in such state C does not depend on Re k ; however, it can still take different values depending on whether turbulence is stationary or decaying.In fact, in the latter case the observed C was ca.twice higher than in the forced case.This observation concerned final stages of decay in the grid turbulence, where the self-similar state was reached.Another conclusion follows for the region closer to the grid, where C Þ const and turbulence is out of equilibrium.
According to Rubinstein and Clark (2017) such a state occurs when turbulence evolves rapidly from one self-similar state to another.Those authors propose to express the energy spectrum as a sum of equilibrium and nonequilibrium parts with Ẽ , , E 0 and where the (possibly nonstationary) E 0 (k, t) satisfies the equilibrium relation Eq. ( 5).Bos and Rubinstein (2017) derived the Kolmogorov's relation for E 0 in the inertial range where C K ≈ 1.5 is the Kolmogorov's constant, and the following form of the nonequilibrium part where is the time derivative of the dissipation rate.Equations ( 7) and ( 8) are assumed to be valid for a wavenumber range k L , k , k h , see Fig. 1.Outside this range E(k, t) was assumed to be zero.Explicit form of E(k, t) was also derived in Yang et al. (2018), using somewhat different approach, which includes dissipative part of the spectrum.Derivation of the nonequilibrium dissipation scaling laws was further performed in Bos and Rubinstein (2017) by assuming that all the considered turbulence properties can be, analogously, decomposed into the equilibrium and nonequilibrium parts, that is, 5 0 1 , k 5 k 0 1 k, L 5 L 0 1 L, and C 5 C 0 1 C and calculated based on the equilibrium and nonequilibrium spectra, Eqs. ( 7) and (8), respectively.This, among other results, leads to the expressions (see appendix A for details) where the last inequality implies that 5 0 1 ≈ 0 .
Using the Taylor series expansion, with k/k 0 being a "small" argument, those authors arrived at the following, approximate result for C : where Re k0 5 k 0 U 0 /n is the "equilibrium" Taylor-scale-based Reynolds number.Equation ( 11) is very close to the experimental result in (3), as the power coefficients n and m in Eq. ( 3) are close to 1. Similarly, the ratio L/k derived by Bos and Rubinstein (2017) reads The power coefficient 1/14 is small and within the accuracy of experimental results, Eq. ( 4) is valid.
We note here that in Bos and Rubinstein (2017) the case of decaying turbulence is studied, as such setting is usually investigated in laboratory experiments.However, those considerations remain valid for another form of nonequilibrium, where the production is locally large and turbulence develops.In such case / in the Eqs.( 8) and ( 9) is positive and so is the nonequilibrium part of the spectrum and k.In atmospheric turbulence both / , 0 and / .0 can be observed.In the former case, Eq. ( 9) predicts (1 1 k/k 0 ) , 1 and it follows from Eq. ( 11) that When turbulence develops, we have (1 1 k/k 0 ) . 1 and obtain C , C 0 : All considerations of Bos and Rubinstein (2017) were performed under the assumptions of homogeneity and local isotropy.These assumptions are usually employed to calculate turbulence dissipation rate from in situ measurements performed by aircrafts, where it is also assumed that the energy spectrum follows the 25/3 law.We also use the assumption of local isotropy in our analyses; however, following Bos and Rubinstein (2017) we account for possible deviations from the Kolmogorov's scaling due to nonstationarity.The next step, left for further work, would be to account for inhomogeneity and anisotropy.This could be done along the line of the studies of Chen and Vassilicos (2022), where transport equations for the second-order structure functions are considered.We note that the nonequilibrium scaling relations (11) and ( 12) are also present in the nonhomogeneous flows, like, e.g., turbulent boundary layers or wakes, cf.Obligado et al. (2016Obligado et al. ( , 2022)), Nedić et al. (2017); hence, we could assume that the dependence of C and L/k on Re k will not be much affected by anisotropy in the atmospheric boundary layer flows, although the proportionality constants may be somewhat different.

b. Estimation of C from turbulence models
Bos (2020) proposed a rough estimate of the upper bound of C , based on common turbulence models.For this, Eq. ( 9) was used: where was replaced by 0 , as the nonequilibrium correction / is small.The term can be estimated from the k-closure equations [cf.Pope (2000) and appendix B], which model the time evolution of the kinetic energy and .It is to note, however, that the transport equation for is purely empirical and its form cannot be derived from first principles; therefore, it should be used only due to a lack of better options.
We extend here the analysis by Bos (2020) by including the buoyancy into the considerations.We assume the possibly simplest form of k-model for buoyant flows, as suggested by Rodi (1987) and assume that transport terms of are negligible in comparison to other components in the equation.With this we have 0 where P is the shear production of the turbulence kinetic energy, G is the buoyancy production (see appendix B), and the model constants take the following values: A 1 5 1.44,A 2 5 1.92, and A 3 5 1.14.It should be noted that the latter appears to be highly case dependent and different values and modifications of the model were reported in the literature.The one used here was suggested by Baum and Caponi (1992).Substituting ( 14) into ( 13) and ( 11) we obtain We will now introduce the flux Richardson number Ri f 5 2G P , which is negative for the unstable and positive for the stable stratification, such that Eq. ( 15) now reads There has been a lot of discussions in the literature about the existence of a certain critical Richardson number above which turbulence is suppressed (Galperin et al. 2007;Zilitinkevich and Baklanov 2002).We refer here to the work of Grachev et al. (2013), who report that the Kolmogorov's 25/3 law fails when Ri f exceeds a critical value 0.2-0.25.However, some small-scale "supercritical" turbulent motions can still survive up to Ri f 5 0.5.We can now consider different values of P/ 0 and Ri f and estimate corresponding C /C 0 from Eq. ( 16).Ri f 5 0.5 will be considered as the upper limit of the flux Richardson number.As reported by Bos (2020) for the case G 5 0, when the production is equal to the dissipation P 5 0 C ≈ 1:1C 0 , which seems to be a good estimate, as the k-model considered predicts stationary state for production slightly exceeding dissipation.For the case without buoyancy, C is bounded from above by the value (Bos 2020) C ≈ 1:8C 0 , which correspond to the case with no shear production P/ 5 0. This follows exactly from Eq. ( 15), with G/ 5 0. We estimated the ratio C /C 0 for nonzero buoyancy from Eq. ( 16).Results are presented in Fig. 2.
Generally when turbulence becomes locally stronger, C is smaller than its equilibrium counterpart.The opposite is true when turbulence becomes locally weaker.It is to note that the value C 5 1.8C 0 remains an upper bound also in the presence of negative buoyancy flux (positive Ri f ), which follows from the fact that (A 1 2 A 3 Ri f ) in Eq. ( 16) remains positive for the maximal value of Ri f 5 0.5.It could, however, be possible that some limited regions of zero turbulence production Values of C /C 0 estimated from Eq. ( 16).
P and negative buoyancy G , 0 would be present.We can expect turbulence would be in a strong nonequilibrium, decaying fast in such zones.Larger values of C were obtained, e.g., in the study of Nedić et al. (2017) from direct numerical simulations data of boundary layer flows at relatively low Reynolds number.However, the derivations of Bos and Rubinstein (2017) were performed for simplified high-Re form of the spectrum, cf.Fig. 1.Hence, the results may not be valid at low Re, where the dissipative part of the spectrum is nonnegligible.Low values of Re k are unlikely to be found in the atmospheric boundary layer turbulence; hence, we would rather expect C to remain bounded from above by 1.8C 0 .

Data and methods
To verify, whether the presented above analysis can be applied to atmospheric turbulence, we used high-resolution data collected in July 2017 during the ACORES campaign in the eastern North Atlantic around the island of Graciosa (Siebert et al. 2021).Measurements were performed with the Airborne Cloud Turbulence Observation System (ACTOS) (Siebert et al. 2006) mounted 170 m below the BO-105 helicopter.We analyze the vertical wind velocity w corrected for platform motion and attitude (cf.Nowak et al. 2021) provided by the ultrasonic anemometer-thermometer (Siebert and Muschinski 2001).The sensor was mounted on ACTOS and has the sampling frequency of 100 Hz.The standard deviations due to uncorrelated noise for velocity measurements are 0.002 m s 21 .Relatively low velocity of the helicopter (∼20 m s 21 ) and high frequency of the sensor allows us to reliably resolve small scales of turbulence down to ∼0.5 m.The further analysis is based on the vertical velocity component only, as the horizontal components still showed some influence of platform attitude and motion which could not be completely removed with standard procedures.These issues caused problems with Reynolds decomposition and subsequent analysis.We note that Obligado et al. (2022) compared estimates of C with U based on turbulence kinetic energy with estimates where only one component of velocity was used to calculate U.In spite of a slight shift of the data points due to anisotropy of the flow in the turbulent boundary layer, the scaling remained unaffected.We will use here the latter methodology and define U 5 w 2 : Helicopter flights during ACORES were performed over the ocean inside the 10 km 3 10 km square adjacent to the Graciosa island.The typical flight duration was 2 h.The trajectory included vertical profiles up to 2000 m and horizontal legs of the length of several kilometers.In our analysis we use data from the horizontal segments for the same flights as selected in Nowak et al. (2021): flight 5 on 8 July 2017 and flight 14 on 18 July 2017, distinctive by the stratocumulus presence and STBL stratification (considerably well mixed in flight 5, considerably decoupled in flight 14).Detailed analysis of turbulence statistics based on these data was performed in Nowak et al. (2021).The study of nondimensional dissipation coefficient and L/k ratio performed in this work provides, additionally, estimation of turbulence temporal tendencies in the two types of STBL.

a. Modified spectra
Due to finite frequency of ACTOS turbulence sensor, the wind velocity is measured with a spectral cutoff.As mentioned, ACORES data have the spatial resolution of ∼0.5 m, which is still larger than the typical size of the smallest Kolmogorov's eddies (∼0.001 m to ∼0.01 m).With this, part of the energy spectrum remains unresolved and to recover value of indirect methods, which rely on the Kolmogorov's hypotheses, should be used.
In this subsection we show that in spite of the presence of the nonequilibrium correction, the turbulence kinetic energy dissipation rate can still be calculated from the onedimensional spectra by fitting to the 25/3 slope, provided that the fitting range is moved toward larger wavenumbers.This is possible, because the nonequilibrium correction affects the one-dimensional spectra only slightly and its influence is negligible for large wavenumbers.Hence, we can still estimate from the measured velocity signals even in the presence of the spectral cutoff at the Nyquist frequency f s /2 5 50 Hz.
In the isotropic turbulence, relation between onedimensional longitudinal spectra and the spectral energy density reads (Pope 2000) E(k, t) is described in Eq. ( 6) and its equilibrium and nonequilibrium parts are given in Eqs. ( 7) and ( 8): The ratio / ≈ / 0 can be estimated from Eq. ( 14) for selected values of P/ and G/ and the energy spectrum can be determined if k 0 and 0 are known.Plots of the onedimensional spectra E 11 , calculated numerically from Eq. ( 17) after substituting Eq. ( 18), corresponding to the case where P/ 5 2, G/ 5 0 (developing turbulence with the production twice larger than the dissipation) and P/ 5 0, G/ 5 0 (decaying turbulence, no turbulence production) are presented in Fig. 3.We used k 0 5 0.06 m 2 s 22 and 0 5 5 3 10 24 m 2 s 23 , which are sample values (to the order of magnitude) measured in the STBL.
To estimate TKE dissipation rate, we assume that within a certain range of wavenumbers the energy spectrum follows the 25/3 law where C 11 ≈ 0.49.Linear least squares fit procedure in the selected range of wavenumbers is performed to estimate the coefficient and calculate .The above formula is no longer true in the case of modified spectra; however, if the fitting range is moved toward larger wavenumbers, the influence of Ẽ should be relatively small, because this term decreases faster with k than E 0 .We used here the range of wavenumbers [1 2 10] m 21 .Our reference case will be P 5 , G 5 0, for which we assume that E 11 follows Eq. ( 19).For this case the linear least squares fit procedure provides ref 5 4.55 3 10 24 m 2 s 23 , which is somewhat smaller than the value prescribed a priori (equal to 0 5 5 3 10 24 m 2 s 23 ).This follows from the fact that Eq. ( 14) predicts / somewhat smaller than zero for P 5 .We next modified P/ and G/ and calculated the corresponding in the same range.Results of / ref and C /C 0 determined from Eq. ( 15) are given in Table 1.As it is seen, even in the case of production being locally 3 times larger than the dissipation, the overprediction of is around 10%.Analogously, in case of no production is underpredicted by 7%.Most of the values of C /C 0 , calculated for ACORES data, are within the range 0.6C 0 -1.8C 0 (see section 4).Hence, as it follows from the above estimation, we expect the over/underprediction of due to modification of the spectra should be less than 10% for the chosen fitting range.To minimize this error it is desirable to move the fitting range toward possibly largest resolved wavenumbers, but such that the effect of aliasing due to finite sampling frequency is still negligible.Wacławczyk et al. (2020) argued that estimates based on the structure function approach are better in such case.This is also confirmed in the experience of Siebert et al. (2006) and Nowak et al. (2021).For this reason, in the subsequent analysis we estimate using the transverse second-order structure function approach, by fitting the data to the formula D 2 ≈ C 2 (r) 2/3 , where C 2 5 2.8. Figure 4 presents structure functions estimated from airborne measurements in the decoupled STBL.The wind velocity was recorded at height 230 m in a region of large turbulence production and at 990 m (in a region inside the cloud with locally weak turbulence production).As it is seen S 2 profiles clearly deviate from the Kolmogorov's scaling at high r.

b. Characteristic time scales
For each configuration it is also possible to estimate the characteristic time scale, which we define as The term is present directly in Eq. ( 18) for the modified form of the energy spectrum; hence, we can expect that 1/t defined above is a measure of nonequilibrium.The equilibrium time scale of turbulence will be defined by T 0 5 k 0 / 0 ≈ k 0 /.
To estimate the ratio t/T 0 we use Eq. ( 14) which is rewritten below as and use P/ and G/ as prescribed in Table 1.Results are presented in Table 1, last column.As it is seen, the further the system is from its stationary state, the smaller the time scale t becomes.Hence, t can be interpreted as the time scale of turbulence "adjustment" toward equilibrium.

c. Effects of Reynolds decomposition
To calculate turbulence statistics the instantaneous velocity should first be decomposed into the mean and fluctuating 3. One-dimensional spectra calculated for the case P/ 5 2, G/ 5 0 (developing turbulence with the production twice larger than the dissipation): solid red line; P/ 5 0, G/ 5 0 (decaying turbulence, no turbulence production): dot-dashed blue line; compared with K41 spectrum: dashed black line.parts.Here, the ensemble average is approximated by the space average along the flight track.The size of averaging window used for detrending, AW D , should be much larger than the characteristic turbulence length scale but much smaller than the length associated with the mean changes.In practice, in the atmospheric turbulence, the presence of largescale convective motions, internal waves, and changes of atmospheric conditions along the flight track makes the choice of the averaging window difficult.In this work we perform the analysis for vertical velocity components, for which the changes of the mean component were of lesser magnitude.In Nowak et al. (2021), AW D 5 50 s, which correspond to, approximately 1 km length, were used for the horizontal segments.Such window was approximately 10 times larger than the estimated integral turbulence length scales which were of order 100 m.
Size of the averaging window for detrending AW D affects mostly the estimation of the integral length scale and the turbulence kinetic energy.We investigated how the values of dimensionless dissipation C determined from the formula are sensitive to the choice of AW D .The integral length scale L 5 L 11 was calculated from the transverse two-point correlation coefficient as described in appendix C, was estimated from the second-order structure function in the fitting range [0.6-6] m, which corresponds to the range of wavenumbers considered in section 3a.Finally U was calculated as the standard deviation of the vertical component of fluctuating velocity w.
The analysis is performed for the subcloud segment F05LEG307 from flight 5, characterized by low values of C , which suggests nonequilibrium, developing turbulence, a part of subcloud segment F14LEG287 from flight 14, where turbulence is close to equilibrium and cloud segment F14LEG992 from flight 14 with C .C 0 .For the first two segments the conditions remain approximately constant along the flight track.In the third segment, larger variations of turbulence properties along the flight track are observed; however, the mean value of C still exceeds the equilibrium value, see section 6.After detrending another averaging window AW S should be chosen for the calculation of statistics.This window should be sufficiently long to reduce the bias and the random errors to acceptable levels (Lenschow et al. 1994).In the following subsection, we estimate the size of AW S sufficient to estimate C with a good accuracy.Here, we perform averaging over the whole considered parts of the segments, that is, AW S 5 270 s for flight 5 and AW S 5 300 s for both segments of flight 14.These lengths correspond to at least ≈ 30L for the largest calculated L.
Results are presented in Table 2.As can be seen, the larger AW D is chosen, the larger resulting L and U are.By increasing AW D even further we would get slower growth and eventually convergence of L and U.In spite of the dependence of L and U on AW D , changes of C are much smaller.E.g., once AW D is reduced from 75 to 50 s, C changes by 0.01 or 0.02.If AW D is further reduced, the value of C becomes close to equilibrium.This would be expected, as both L and U are influenced by large scales, which are filtered out if AW D is too much reduced.We decided to use AW D 5 50 s, the same as in Nowak et al. (2021), for further analysis.Such value still allows us to detect differences in C levels, and, with moderate L, the size of averaging windows for the calculation of statistics, AW S , could be taken small enough to observe variations of turbulence properties along flight tracks.

d. Effect of averaging windows
Keeping AW D 5 50s constant we next investigated influence of the results to the choice of AW S , that is, the window used to calculate statistics from the previously detrended signal.The analysis is performed for the first half of F14LEG287 from flight 14 in the decoupled STBL, where C remains approximately constant.The largest AW S 5 200 s corresponds to, approximately 32L, or ≈4 km.We moved the window by 5 s (≈100 m) along the flight track and calculated turbulence statistics.Such procedure allows us to detect possible variations of atmospheric conditions along the flight track.For this, it would be optimal to have possibly small AW S , as otherwise changes of C are averaged out.On the other hand, AW S should be large enough to calculate statistics with a good accuracy.We gradually decreased AW S to 150, 100, and 50 s which corresponds to, approximately 25L, 16L, and 8L, respectively.Results are presented in Fig. 5. Standard deviations std(C ) for the investigated segment increase with the decreasing AW S and equal, respectively, std(C ) 5 0.036, 0.045, 0.064, and 0.21.The last number is clearly unacceptable (which is also visible in Fig. 5), as it should be possible to detect changes of C on the order of 0.1 to draw conclusions about the state of turbulence in the investigated region.We decided to further use such averaging window AW S which correspond to the length larger than 20L.

Results: Analysis of turbulence states in STBL
We first address and compare two horizontal flight segments of comparable, low altitudes of 287 m (decoupled STBL, F14LEG287) and 307 m (coupled STBL, F05LEG307).In the former case the averaging window AW S 5 150 s was used, in the latter, the window was increased to AW S 5 170 s due to larger length scales.The window was moved every 5 s (≈100 m) along the flight track and each time turbulence statistics were calculated.Figure 6a presents values of C as a function of the position of the center of the averaging window.In the figure two parallel lines of a constant C 0 5 0.45 and C 5 0.8 are plotted.The first one would correspond to the equilibrium, stationary case.We note that this equilibrium value may in fact differ slightly for different flows, as it is weakly influenced by the large turbulence structures and their scaling (Valente and Vassilicos 2012;Bos et al. 2007).This influence was not taken into account in the theoretical considerations in section 2, as it was assumed that the energy density has sharp spectral cutoffs at k L and k h .Moreover, possible flow anisotropy may also play a role and modify the proportionality constant.We choose C 0 5 0.45 as it is the value obtained for F14LEG287 at small AW D 5 15 s, see Table 2, for which nonequilibrium effects were mostly filtered out, as it was discussed in section 5d.The second value, C 5 0.8 corresponds to a freely decaying turbulence with P 5 0, where the equilibrium, self-similar state was reached (see Table 1).
As it is seen in Fig. 6a, initially, C ≈ C 0 .As the averaging window is moved further along the flight track, a moderate increase of C is first observed, followed by a sharp decrease toward smaller values.It would suggest that a weak turbulence decay, followed by a region of strong turbulence production was detected.In the analogous figure for the coupled STBL, Fig. 7a, only C , C 0 are detected, which can be again interpreted as a region of strong turbulence production.
We next compared two subcloud flight segments: F14LEG448 in the decoupled STBL of height 448 m and F05LEG553 in the coupled STBL of height 553 m.In Nowak et al. (2021) nearly zero buoyancy production at those levels was reported.The weak turbulence production is manifested through larger average values of C , both in decoupled and coupled case, cf.Figs.6b and 7b, respectively.Insufficient turbulence production at this altitude may lead to decoupling, when the eddies fail to mix air over the entire depth of the STBL, and consequently STBL separates into two parts: cloud driven at the top and surface driven at the bottom, cf.Nowak et al. (2021).
Finally, we consider in-cloud segments F14LEG992 in the decoupled STBL and F05LEG819 in the coupled STBL.The signals were recorded at heights 992 and 819 m.In both cases values of C vary considerably, cf.Figs.6c and 7c; however, they are somewhat larger in the decoupled case which indicates weaker turbulence production and is again in line with the observations of Nowak et al. (2021).The results of C can be compared with experimental estimations of Zheng et al. (2021) in high-Re turbulence behind a grid, where C first decreased in the production region and next increased in the nonequilibrium decay region toward value C ≈ 0.8.We also mention here the direct numerical simulation study of Gallana et al. (2022), who calculated C in stably and unstably stratified flows.In the former case C increased toward C ≈ 0.8, indicating turbulence suppression and in the unstable conditions decreased toward C ≈ 0.4, below the equilibrium value.In STBL turbulence is manly driven by radiative and evaporative cooling at the cloud top, leading to instability.However, boundary layer is also capped from above by a stably stratified layer.This may explain large variations of C along the in-cloud flight tracks.
To classify whether turbulence is in its equilibrium or nonequilibrium state, we calculated the length scales and the local Reynolds numbers and plotted L/k and C as a function of Re k .The data points were compared with the theoretical predictions (lines).For C , in the equilibrium case in high-Re atmospheric flows, we should obtain C 5 const.In the nonequilibrium, dependence of C on Re k should be described by Eq. ( 11).As we do not have the knowledge about the equilibrium value Re k0 , the coefficient in Eq. ( 11) was calculated from the linear least squares algorithm.In the lowest flight segment F14LEG287 in decoupled STBL statistics calculated at x , 4500 m follow, approximately the equilibrium line C 5 const (black dots in Fig. 8a).However, when the position of the center of the averaging window exceeds 4500 m the nonequilibrium states can be detected (red stars in Fig. 8a).The values of C calculated for coupled STBL for F05LEG307 decrease slightly with Re k , which suggests the presence of nonequilibrium scaling (cf.Fig. 9a).
Figures 8b and 9b refer to the flight segments F14LEG448 and F05LEG553 below the cloud.Figure 8b suggests that part of the signal in the decoupled STBL is a record of equilibrium turbulence, while statistics of another part fit the nonequilibrium formulas.For the coupled STBL the statistics follow rather the equilibrium predictions, cf.Fig. 9b.A large scatter of results is visible in Fig. 9b such that it is not easy to draw conclusions.Note that the flight segment F05LEG553 is relatively short and so is the range of the detected Re k values.Equilibrium scaling may suggest that in spite of the weak turbulence production, the changes of atmospheric conditions are relatively slow, such that the turbulence spectra already relaxed to the self-similar form (5). Similar, nearly constant levels of C were reported, e.g., by Neunaber et al. (2022) in the study of turbulence in wakes behind turbines.
Particularly interesting is the study of the in-cloud segments F14LEG992 and F05LEG819.In the decoupled STBL most of the points follow the nonequilibrium scaling relations, cf.We may call such state of the system as a quasi-equilibrium state, which means that the changes of external conditions are slow enough, such that the system undergoes through several local equilibria.
We note here that for certain points the assignment to equilibrium or nonequilibrium region may be somewhat arbitrary.This concerns especially those points which represent statistics calculated at the boundary between the equilibrium and nonequilibrium region.Sometimes this assignment becomes more definite when the dependence of the length scale ratio L/k on Re k is studied, see Figs. 10 and 11.In the equilibrium case, the statistics should follow Eq. ( 2), while in the nonequilibrium L/k and Re k are inversely proportional and follow Eq. ( 12).Both functional dependencies are plotted in Figs. 10 and 11 as the black and red lines, respectively.We note that we did not perform another curve fitting here, the constants are those estimated for data from Figs. 8 and 9 for respective flight segments and divided by 15, as follows from the Eqs.( 11) and ( 12).The good agreement confirms validity of the theoretical derivations of Bos and Rubinstein (2017).After the classification of statistics in Figs.8-11 the nonequilibrium parts were marked in Figs. 6 and 7 (as red rectangles).
Nonequilibrium scaling relations were found in many laboratory studies, as examples we mention here the works of Zheng et al. (2021), which concerns high-Re flows behind a grid, and of Nedić et al. (2017), which concerns boundary layer flows.
We next considered separately F14LEG143 of the decoupled STBL of altitude 143 m, as there are no measurement data of comparable height for the coupled STBL.As it is seen in Fig. 12, two region of C Þ C 0 can be distinguished.The part on the left, where C .C 0 suggests turbulence decay.Statistics calculated at the final part of the signal, from F14LEG143, with C , C 0 , imply on the other hand the dominant role of the production.On the middle part of the curve in Fig. 12, C oscillates around the equilibrium value of 0.45.These equilibrium values are marked as black dots in Fig. 13a.Apparently, they do not depend on Re k .Moreover, the ratio L/k increases with Re k , as predicted by the equilibrium relation (2), see Fig. 13b.However, the remaining points, marked as red and blue stars in Fig. 13, do not follow one nonequilibrium line.We can argue that they correspond to two distinct regions, with different initial conditions and different Re k0 .By using the least squares algorithm we calculated two coefficients of the nonequilibrium relation ( 11 (blue stars in Fig. 13).We found Re k0B ≈ 1.2Re k0A .Assuming C 0 ≈ 0.45 it is possible to recover the initial-state Re k0 from the calculated A and B coefficients.We found Re k0A ≈ 7500 and Re k0B ≈ 9000 as the equilibrium (initial) Reynolds numbers in the respective regions.We next plotted the red points for rescaled coordinate Re k 5 1:2 Re k and found, as expected, that they follow the same nonequilibrium curve as the blue points, see Fig.

Conclusions
The main purpose of this work was to show that certain information about the temporal variations of turbulence and the budget of the turbulence kinetic energy in the STBL can be acquired by analyzing statistics of wind velocity based on in situ measurement data.The key indicators are the nondimensional dissipation coefficient C and the integral-to-Taylor-scale ratio L/k.C takes largest values, up to around C ≈ 0.8 when turbulence decays, reaches its equilibrium C ≈ C 0 ≈ 0.45 in the stationary state and becomes smaller C , C 0 when turbulence becomes stronger.We showed here that the upper limit C ≈ 1.8C 0 should not be exceeded even under a stable stratification, due to the presence of the critical Richardson number, above which turbulence is suppressed.
Moreover, as it was discussed in the previous works (Bos and Rubinstein 2017) two different states of turbulence system can be identified: equilibrium, where the turbulence kinetic energy spectrum has a self-similar shape, (5), and a nonequilibrium where a correction to the spectrum should be considered; see Eq. ( 6).The equilibrium state can possibly be nonstationary; however, the changes of external conditions (forcing) are slow enough such that the system has enough time to relax and the energy spectrum takes the self-similar form (5). To detect the nonequilibrium states we studied dependence of C and L/k on the Taylor-based Reynolds .This method complements measures of nonstationarity and nonequilibrium reported by Mahrt and Bou-Zeid (2020, see also references therein), who mention the lack of an inertial subrange and results of direct production-dissipation balance estimates.The latter are not available from horizontal flight segments, as derivatives in the vertical direction cannot be calculated.However, as it was discussed in section 2b, the production-dissipation balance is related to the ratio C /C 0 .
In this work we focus on the analysis of airborne measurement data, which often suffer from various deficiencies, related to finite frequencies of the sensors or short available averaging windows (Wacławczyk et al. 2017;Akinlabi et al. 2019).We showed that despite these deficiencies both C and L/k can still be estimated from the ACORES data (Nowak et al. 2021).In particular, to estimate it is crucial to move the fitting range toward smaller scales (larger wavenumbers) to minimize the influence of nonequilibrium correction to the spectrum.This 27/3 correction is generally subdominant and influences larger scales.We also showed that the indicator functions can be estimated with a sufficient accuracy if the size of the averaging window used to calculate statistics is larger than 15L.This allows us to obtain C and L/k for horizontal flight segments.Results compare very well with the theoretical predictions of Bos and Rubinstein, which concerns both dependence on local Re k as well as on the initial conditions Re k0 , which is shown by a proper rescaling of the data.Moreover, with a curve fitting of the data we were able to recover values of Re k0 .
The analysis performed here allows for additional comparison of turbulence in the coupled and decoupled STBL's, previously studied by Nowak et al. (2021).In both cases, C is smaller at lower altitudes due to strong shear and/or buoyancy production.On the other hand, in the subcloud parts C increases above C 0 which suggests that in these regions turbulence decays.Large variations of C were detected inside the cloud in both cases.Apparently, along the flight track, regions where turbulence production surpass or is in balance with the dissipation alternate with areas of turbulence decay.Here, the main difference between the coupled and decoupled STBL was observed.Namely, in the former, turbulence seems to undergo a series of equilibrium states, where C ≈ const and L/k ∼ Re k , while in the latter case a clear nonequilibrium is observed.This may be a consequence of a different flow organization in the latter case, with smaller and faster large-scale circulations (Nowak et al. 2021).Moreover, at all the investigated altitudes, mean C was higher in the decoupled than in the coupled case which indicates the weaker turbulence production in the decoupled STBL.
In this work only limited amount of data from two ACTOS flights were investigated.In future, more thorough analysis to better characterize turbulence in the coupled and decoupled STBL is foreseen.The same analysis can be repeated for various measurement data in the atmospheric boundary layers where the impact of transient states is nonnegligible (cf.Mahrt and Bou-Zeid 2020), e.g., to study the collapse of turbulent convective daytime boundary layer (El Guernaoui et al. 2019;Wang et al. 2020;Lothon et al. 2014), in and between various types of the clouds in order to better understand their dynamics and even in clear air turbulence to characterize its dynamic properties.
Another interesting aspect is the relation of the present study to the K62 theory (Kolmogorov 1962), which takes into account temporal and spatial variations of the instantaneous dissipation rate.This phenomenon is known as the "internal intermittency" and affects the scaling of higher-order structure functions.A careful analysis of the fine-scale structure of turbulence performed by Siebert et al. (2010) based on ACTOS data revealed the intermittency coefficient in cloud turbulence is consistent with values obtained in laboratory experiments.A question to be asked is how the nonequilibrium structure functions are affected by the small-scale intermittency and whether temporal variations of statistics considered here can be linked to the premises of the K62 theory.
The presence of nonequilibrium in the stratocumulus clouds indicate that the common turbulence closures like the Smagorinsky LES model may fail to predict the dynamics of such clouds correctly.In particular, this may influence the prediction of boundary layer decoupling or cloud dissipation.A nonequilibrium turbulence model may largely improve these predictions.This is another direction of future studies.parts, that is, 5 0 1 , k 5 k 0 1 k, L 5 L 0 1 L and calculated with the use of Eqs.(A1) and (A2) as follows: and where I 0 5 k 21 E 0 (k, t)dk.Further, k, , and L are calculated analogously, by just replacing E 0 by Ẽ in Eqs.(A3)-(A5).This leads to the following, approximate results: ) and ) where the last inequality follows from the fact that k L , , k h .Further, as L 5 3pI/4k, with the use of Eq. (A5) we can express L as With this, Bos and Rubinstein (2017) proposed to write C and Re k in the following form: where it was assumed that ≈ 0 due to the inequality / 0 , , k/k 0 ; see Eq. (A11).Re k0 above denotes the equilibrium Reynolds number.Next, after substituting (A10) and assuming that Ĩ/I 0 is a first-order term of the Taylorseries expansion, the final formula for the nonequilibrium C was obtained as follows: Similarly, the analog of Eq. ( 4) was presented in Bos and Rubinstein (2017).To derive it, let us note that is the buoyancy flux.Here, b stands for the fluctuations of the buoyancy.The vertical component of the fluctuating velocity u 3 will also be denoted by w; hence, the alternative notation for the buoyancy flux is G 5 wb .
In the two-equations k-closures a gradient diffusion hypothesis is used to model the turbulent transport and the production terms, while another equation is solved for the turbulence kinetic energy dissipation rate (Jones and Launder 1972).This equation, adopted in Rodi (1987) for the study of buoyant flows reads The model constants C m 5 0.09, A 1 5 1.44,A 2 5 1.92, and s 5 1 are standard, the value of A 3 5 1.14 used in this work was proposed by Baum and Caponi (1992).

APPENDIX C
Estimation of the Integral Length Scale from Measurement Data The longitudinal and transverse two-point velocity correlations coefficients are defined as (Pope 2000) f (r, t) 5 1 u 2 L u L (x, t)u L (x 1 e L r, t) (C1) and g(r, t) 5 1 u 2 N u N (x, t)u N (x 1 e L r, t) , (C2) where the subscript L denotes longitudinal vector component (i.e., along the flight track) and the subscript N the transverse component (perpendicular to the flight track).In the isotropic turbulence u 2 L 5 u 2 N and the following relation between f(r, t) and g(r, t) holds (Pope 2000): The longitudinal integral length scale is defined as L 11 (t) 5 f (r, t)dr: (C4) In the isotropic turbulence L 11 ≈ L, where L was applied in Eq. ( 1) and is defined with the use of the spectral energy density function We make use of this formula to estimate L from the measurement data.For this the transverse two-point correlation coefficient g(r, t) is calculated by averaging product of velocity components recorded at two points placed at a distance r from each other.We gradually increase r and obtain an approximate form of the function g. Results are deteriorated for large r, due to small size of the sample; however, we expect g(r, t) is calculated with a good accuracy for small r.Next we find the distance r 5 r 0 at which g(r, t) crosses 0 for the first time.We integrate g(r, t) from 0 to r 0 .The calculated value is next divided by 0.57 and the result should be an estimation of the integral length scale L as predicted by Eq. (C8).

FIG. 1 .
FIG.1.Model for spectral energy density in the equilibrium (black dashed line) and in nonequilibrium (solid red line) considered inBos and Rubinstein (2017).
FIG. 4. One-dimensional transverse structure functions estimated from airborne measurements, from the vertical wind velocity component recorded at height 280 m (region of large turbulence production): solid red line; at height 990 m (region inside the cloud with locally weak turbulence production): dot-dashed blue line; compared with K41 spectrum: dashed black line.Vertical lines mark the fitting range.
FIG. 5. C /C 0 calculated for different AW S from part of signal from flight 14, F14LEG287.Horizontal coordinate x denotes the position of the center of the averaging window.
14a.The same can be done for the relation L/k; however, in this case L/k ∼ Re 21/14 k and to obtain the same coefficient we have to plot the data for differently rescaled Re k 5 (1:2) 15 Re k .Results are plotted in Fig. 14.They again confirm the validity of the analytically derived formulas in (11) and (12).

FIG. 14 .
FIG. 14.As in Fig. 13, but for rescaled Re k .FIG. 13.(a) C as a function of Re k , (b) L/k as a function of Re k .Decoupled STBL, F14LEG143.Solid black lines: equilibrium scaling; dashed red line and dot-dashed blue lines: nonequilibrium scaling, Eq. (11); calculated statistics: symbols.
) again, due to the inequality / 0 , , k/k 0 , the term in the denominator in (A16) is replaced by unity; hence, and (A10), and the use of the Taylorseries expansion, the following relation is obtained: (A14) and the equilibrium relation (2) were used.APPENDIX BTransport Equations for the Turbulence Kinetic Energy and Dissipation RateExact form of the transport equation for the turbulence kinetic energy can be derived from the Navier-Stokes equations and reads(Pope 2000) ) where the turbulent viscosity n T readsn T 5 C m k 2 :

TABLE 1 .
Parameters estimate for different values of P/ and G/.

TABLE 2 .
Estimates of turbulence statistics for different AW D .