Investigation of non-equilibrium turbulence decay in the atmospheric boundary layer using Doppler lidar measurements

This work concerns analysis of turbulence in the atmospheric boundary layer (ABL) shortly before and after sunset. Based on a large set of Doppler lidar measurements at rural and urban sites, we analyze frequency spectra of vertical wind at different heights and show that they increasingly deviate from Kolmogorov's $-\mathrm{5}/\mathrm{3}$ prediction in the measured low-wavenumber part of the inertial range. We find that before sunset, the integral length scales tend to decrease with time. These findings contrast with a classical model of equilibrium decay of isotropic turbulence, which predicts that the scaling exponent should remain constant and equal to $-\mathrm{5}/\mathrm{3}$ and the integral length scale should increase in time. We explain the observations using recent theories of non-equilibrium turbulence. The presence of non-equilibrium suggests that classical parametrization schemes fail to predict turbulence statistics shortly before sunset. By comparing the classical and the non-equilibrium models, we conclude that the former may underestimate the dissipation rate of turbulence kinetic energy in the initial stages of decay.

1 Introduction

Turbulence in the atmospheric boundary layer (ABL) undergoes temporal changes with the diurnal cycle. After dawn and under clear sky, the surface heating produces convection and a boundary layer starts to grow. Shortly before sunset, the convective ABL collapses rapidly, and then a stable nocturnal boundary layer (BL) is formed (Nieuwstadt and Brost, 1986).

During the day, turbulence production due to buoyancy is prevalent. In the afternoon, the buoyancy flux decreases gradually and eventually becomes negative (Sorbjan, 1997). The time when the heat flux crosses the zero level was identified by Nadeau et al. (2011) as the beginning of the evening transition. At this time instant, turbulence in the ABL starts to decay more rapidly than in the afternoon. After sunset, turbulence is still produced by shear and the remaining thermal forcing, albeit mostly in a region close to the surface.

Far enough from the surface, turbulence is usually assumed to be approximately homogeneous and isotropic at scales smaller than the integral length scale, respectively. In spite of this considerable simplification, homogeneous and isotropic turbulence is a subject of ongoing research due to its importance for existing theories (Sagaut and Cambon, 2018). Recently, a number of theoretical works have addressed the parametrization of decaying isotropic turbulence (Vassilicos, 2015; Goto and Vassilicos, 2016; Bos and Rubinstein, 2018). In particular, deviations from the predictions of the Kolmogorov theory have been observed in the initial stages of decay, after the forcing is switched off. Kolmogorov's theory of turbulence is of utmost importance as it is widely used to estimate the turbulence kinetic energy dissipation rate ϵ from measured signals. The dissipation rate determines how fast the kinetic energy of turbulence is transferred into heat at the smallest scales, which are of the order of millimeters in atmospheric turbulence. Such small scales can hardly ever be measured with satisfactory accuracy in the free atmosphere. For this reason, the dissipation rate is estimated indirectly by assuming that the energy injected at large scales by forcing is transported at a constant rate from larger to smaller eddies. This process is known as the energy cascade. Taking this assumption, Taylor (1935) formulated the famous relation between the dissipation rate $\mathit{ϵ}=\mathrm{2}\mathit{\nu }\langle s_{ij}{s}_{ij}{\rangle }^{\mathrm{1}/\mathrm{2}}$, where $s_{ij}=(\partial u_i^{\prime }/\partial x_j+\partial u_j^{\prime }/\partial x_i)/\mathrm{2}$; the turbulence velocity scale $\mathcal{U}={\left(\langle u_i{u}_i^{\prime }\rangle /\mathrm{3}\right)}^{\mathrm{1}/\mathrm{2}}$, where $u_i^{\prime }$ is the ith component of fluctuating velocity; and the characteristic length scale of large eddies ℒ (the integral length scale),

$$\begin{array}{}\text{(1)}& \mathit{ϵ}=C_{\mathit{ϵ}}{\displaystyle \frac{{\mathcal{U}}^{\mathrm{3}}}{\mathcal{L}}},\end{array}$$

where C_ϵ≈0.5. Equation (1) forms the basis of many turbulence parametrization schemes.

The validity of the Taylor law was questioned in recent theoretical works and experimental observations; see Vassilicos (2015). In particular, it was observed that at the onset of decay, the dissipation rate followed a non-standard relation:

$$\begin{array}{}\text{(2)}& \mathit{ϵ}=C_{n\mathit{ϵ}}{\mathcal{U}}_{\mathrm{0}}{\mathcal{L}}_{\mathrm{0}}{\displaystyle \frac{{\mathcal{U}}^{\mathrm{2}}}{{\mathcal{L}}^{\mathrm{2}}}},\end{array}$$

where C_nϵ is a constant and 𝒰₀ and ℒ₀ denote initial values of the turbulence velocity scale and the length scale. Bos and Rubinstein (2018) argued that the appearance of Eq. (2) is connected with deviations from the $-\mathrm{5}/\mathrm{3}$ scaling of the frequency spectra in the low-wavenumber part of the inertial range. These deviations are observed when the energy transfer rate across the scales is not constant due to a sudden change in forcing. This state will be further referred to as “non-equilibrium”. The same notion is used in thermodynamics and fluid mechanics to describe states after a sudden change in external conditions when the system evolves towards another equilibrium (Wacławczyk, 2021). Here, “equilibrium” is related to Kolmogorov's turbulence characterized by the $-\mathrm{5}/\mathrm{3}$ law or C_ϵ≈const. Non-equilibrium, on the other hand, denotes deviations from these laws.

The non-equilibrium scaling has been observed in a number of laboratory and numerical experiments (Valente and Vassilicos, 2012; Obligado and Vassilicos, 2019; Steiros, 2022a; Zheng et al., 2023; Obligado et al., 2022; Steiros, 2022a) as well as in atmospheric turbulence (Wacławczyk et al., 2022). The latter work concerned analysis of airborne measurement data from the stratocumulus-topped boundary layers (STBLs). Non-equilibrium dissipation scaling of a form close to that in Eq. (2) was observed, especially near the surface and inside clouds in the decoupled STBLs. As discussed in (Wacławczyk et al., 2022), weaker turbulence was characterized by larger values of C_ϵ. In particular, C_ϵ tended to be larger in decoupled STBLs, when turbulence was too weak to mix air over the entire height of the BL. This was in line with previous findings by Nowak et al. (2021), who speculated that turbulence in the decoupled STBLs might be decaying.

To the best of the authors' knowledge, the concept of non-equilibrium has not been discussed in the context of the collapse of the convective BL, although previous studies, when analyzed from this new perspective, deliver strong indications of the non-equilibrium turbulence decay. For example, Lothon et al. (2014) analyzed data from the Boundary-Layer Late Afternoon and Sunset Turbulence (BLLAST) field experiment and calculated the integral length scales ℒ_w from the measured vertical velocity. As reported by those authors, ℒ_w initially decreased in the surface layer and next sharply increased after 19:00 local time (LT). In the mixed layer ℒ_w at first remained constant and then started to increase around 17:00 LT. In classical equilibrium turbulence, the integral length scale is expected to increase in time during the decay. In the non-equilibrium decay, on the other hand, it may initially decrease and next increase in time (Steiros, 2022b). This scenario is consistent with results of Lothon et al. (2014) in the surface layer.

Results of the BLLAST experiment were also analyzed by Darbieu et al. (2015) and were compared to results of numerical simulations. In both cases, the authors found deviations from the Kolmogorov $-\mathrm{5}/\mathrm{3}$ law before sunset, although the validity of the Taylor law (Eq. 1) was not discussed therein. El Guernaoui et al. (2019) discussed the time evolution of the kinetic energy spectra in a numerically simulated ABL. They argued that the decay of the kinetic energy is not uniform across the scales and that the largest scales are the most affected.

Dissipation rates in the surface layer during the afternoon and evening transition were reported by Nilsson et al. (2016). The authors investigated a relation similar to that of Eq. (1) but with the dissipation length l_ϵ instead of the integral length scale ℒ. They concluded that assuming l_ϵ∝z, where z denotes the height, is not sufficient to parametrize the dissipation rate in the surface layer. Instead, they proposed to relate l_ϵ to both z and the height of the boundary layer. Lampert et al. (2016) studied the anisotropy of turbulence during the evening transition and reported that the standard deviation of the vertical velocity decreases with time faster than the corresponding standard deviations of horizontal components. This observation could indicate that the non-equilibrium affects the spectra of the vertical velocity component more than the horizontal ones.

The current work focuses on the decay of turbulence before sunset in order to investigate whether it can be parametrized with the non-equilibrium laws. We analyze data from wind Doppler lidar, which has become a strategic instrument in atmospheric research because it provides vertical profiles of the radial wind component with high spatial and temporal resolution. Our aim is to examine deviations from the Kolmogorov scaling during the evening transition at different heights within the ABL and over different environments. Atmospheric turbulence, which is characterized by huge Reynolds numbers, is an ideal test bed for verifying turbulence theories. Additionally, the theories aim to improve turbulence parametrization schemes, which are of importance for numerical weather prediction and climate models.

It is challenging to provide lidar data with sufficient resolution to recognize deviations from Kolmogorov's scaling in the low-wavenumber part of the inertial range. The spectra are also affected by filtering (averaging) in space and the low signal-to-noise ratio. Hence, these effects should be carefully differentiated. To examine deviations from the Kolmogorov scaling, we also analyze the second-order structure functions. They are mathematically equivalent to the one-dimensional spectra but may respond differently to errors due to the finite frequency of measurements and due to spatial averaging (Schröder et al., 2024; Wacławczyk et al., 2017, 2020). In this work we investigate the scaling of both frequency spectra of vertical wind and the structure functions to assess how they change during the decay of turbulence before sunset. Moreover, we calculate the standard deviations of the vertical velocity and integral length scale and study how they change in time. We also compare the dissipation rates predicted by the classical (Eq. 1) and the non-equilibrium relations (Eq. 2).

The diversified large datasets of sufficient resolution investigated in this work offer a unique opportunity to describe the main differences between decay of turbulent flows generated over urban and rural environments. The rural environment resembles quasi-ideal laboratory conditions for a turbulent flow. The atmospheric boundary layer over an urban environment can experience more shear-driven turbulent flows due to the urban friction and its surface roughness and more thermally generated turbulent flows. The latter are caused by the interactions of the surface with solar radiations (Svensson, 2004; Edokpa and Nwagbara, 2018). Of particular importance is the thermal heat capacity of surfaces and the related urban heat island phenomenon, which is an effect of the heat accumulation in and over an urban area (Oke, 1987). As reported by Nadeau et al. (2011), the decay of turbulence kinetic energy scales with the characteristic time of the heat flux decay. This timescale is smaller for surfaces that cool down more rapidly. In our work we compare the results at both sites (rural and urban) to assess how surface heterogeneity and surface heat capacity affect turbulence properties.

The paper is structured as follows. The theory of equilibrium and non-equilibrium decay is presented in Sect. 2. In Sect. 3 experimental sites and instrumentation are described. Section 4 is devoted to methodology; this is followed by Sect. 5, which contains data analysis. Finally, conclusions and perspectives are discussed in Sect. 6.

2 Theory

2.1 Non-equilibrium spectra and structure functions

The theory of Kolmogorov is the foundation of turbulence research (Pope, 2000). It states that at sufficiently large Reynolds numbers and under the assumption of local isotropy, there exist a range of scales where statistics of velocity take a self-similar form. Further, within this range, a sub-range of scales of size r, where $\mathit{\eta }\ll r\ll \mathcal{L}$, can be distinguished. This will be further referred to as the “inertial range”. In this range of scales, statistics of turbulence depend not on viscosity but only on the dissipation rate ϵ. It follows that the wavenumber spectrum of the vertical wind velocity component can be expressed as

$$\begin{array}{}\text{(3)}& E_{⟂}(\mathit{\kappa },t)=C_{⟂}{\mathit{\kappa }}^{-\mathrm{5}/\mathrm{3}}{\mathit{ϵ}}^{\mathrm{2}/\mathrm{3}},\end{array}$$

where $C_{⟂}\approx \mathrm{0.65}$ is a constant and κ is the wavenumber. Equivalently, the same can be presented in terms of the second-order structure function, which for the vertical velocity component w reads

$$\langle \mathit{\delta }{w}^{\mathrm{2}}\rangle =\langle \left(w\right(\mathbf{x}+\mathbf{r},t)-w(\mathbf{x},t){)}^{\mathrm{2}}\rangle .$$

Within the inertial range and under the assumption of local isotropy, this function takes the form

$$\begin{array}{}\text{(4)}& \langle \mathit{\delta }{w}^{\mathrm{2}}\rangle =C{\left(\mathit{ϵ}r\right)}^{\mathrm{2}/\mathrm{3}},\end{array}$$

where C≈2.86. Equations (3) and (4) form the basis of various schemes for the estimation of the energy dissipation rate, including from lidar measurements (O'Connor et al., 2010; Lothon et al., 2009; Sanchez Gomez et al., 2021).

Equations (3) and (4) work well when turbulence is close to isotropic, at least locally, and is stationary. Recently, extensions of Kolmogorov's theory towards unsteady turbulence were put forward by Vassilicos (2015). These extensions predict that the rate of energy transfer across scales in the inertial range is affected by turbulence decay and is not constant. Bos and Rubinstein (2018) expressed the turbulence kinetic energy spectrum as a sum of the equilibrium, Kolmogorov spectrum and a non-equilibrium correction. They derived a formula similar to Eq. (2) and argued that deviations of the spectra from Kolmogorov's scaling are related to the deviations from C_ϵ=const in Eq. (1). Goto and Vassilicos (2016) focused on the large-scale part of the turbulence kinetic energy spectrum during the non-equilibrium decay and found that it has the self-preserving form

$$\begin{array}{}\text{(5)}& E(\mathit{\kappa },t)\propto \mathit{ϵ}{\mathcal{L}}^{\mathrm{3}}f\left(\mathit{\kappa }\mathcal{L}\right).\end{array}$$

Steiros (2022a) introduced the notion of “balanced nonstationary turbulence”, where the transfer across the scales is proportional to dissipation, albeit with a proportionality constant smaller than 1. This assumption led to a modified form of the energy spectrum in the inertial range:

$$\begin{array}{}\text{(6)}& E(\mathit{\kappa },t)=C_k\mathit{ϵ}(t{)}^{\mathrm{2}/\mathrm{3}}{\mathit{\kappa }}^{-\mathrm{5}/\mathrm{3}}{\left[\mathrm{1}-c(\mathit{\kappa }+\mathcal{L}(t){)}^{-\mathrm{2}/\mathrm{3}}\right]}^{\mathrm{2}},\end{array}$$

where c is a dimensionless constant. The spectra followed the above formula even during the equilibrium decay with C_ϵ=const. The function in square brackets in the above equation reaches the value 1 asymptotically, at large wavenumbers (small scales), where the spectra remain close to the Kolmogorov $-\mathrm{5}/\mathrm{3}$ form.

Obligado and Vassilicos (2019) investigated how the inertial range of the structure functions is affected during non-equilibrium decay of turbulence. They concluded that in the case of decaying turbulence, the second- and the third-order structure functions are closest to Kolmogorov's predictions at the small-scale end of the inertial range. For larger scales, the structure functions increasingly deviate from equilibrium, even at very large Reynolds numbers. The authors considered Lundgren's formula for the structure functions, derived with the use of matched asymptotic expansions. For very high Reynolds numbers, it reads

$$\begin{array}{}\text{(7)}& \langle \mathit{\delta }{u}^{\mathrm{2}}\rangle =C{\left(\mathit{ϵ}r\right)}^{\mathrm{2}/\mathrm{3}}\left[\mathrm{1}-A_{\mathrm{2}}(r/\mathcal{L}{)}^{\mathrm{2}/\mathrm{3}}\right],\end{array}$$

where A₂ is a dimensionless constant of the order 1. Under the assumption of local isotropy, the above formula is mathematically equivalent to Eq. (6), provided that the bracketed term in Eq. (6) can be expanded in the Taylor series.

We now discuss turbulence statistics, in particular, the integral length-scale change in time according to the theory of equilibrium and non-equilibrium decay. Based on this, we can identify which type of parametrization more adequately describes the collapse of the convective boundary layer before sunset. Here we mostly refer to the recent papers by Goto and Vassilicos (2016) and Steiros (2022b), who investigated decaying turbulence using numerical experiments and derived non-equilibrium decay laws.

2.2 Equilibrium decay

Under the assumption of horizontal homogeneity, the transport equation for the turbulence kinetic energy $k=\langle u_i^{\prime }{u}_i^{\prime }\rangle /\mathrm{2}=\mathrm{3}/\mathrm{2}{\mathcal{U}}^{\mathrm{2}}$ in the ABL reads (Pope, 2000)

$$\begin{array}{}\text{(8)}& {\displaystyle \frac{\partial k}{\partial t}}=-T-\mathit{ϵ}+P+B,\end{array}$$

where

$$\begin{array}{rl}& T={\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}\left({\displaystyle \frac{\partial \langle p^{\prime }{w}^{\prime }\rangle }{\partial z}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{\partial \langle u_i^{\prime }{u}_i^{\prime }{w}^{\prime }\rangle }{\partial z}}\right),\\ & P=-\langle u_i^{\prime }{w}^{\prime }\rangle {\displaystyle \frac{\partial \langle u_i\rangle }{\partial z}},\\ & B=\langle w^{\prime }{b}^{\prime }\rangle .\end{array}$$

T, P and B stand for the turbulent transport, shear production and buoyancy forcing, respectively. Above, p^′ denotes the pressure fluctuations and b^′ the fluctuations in buoyancy. The buoyancy flux will be defined as

$$\begin{array}{}\text{(9)}& \langle w^{\prime }{b}^{\prime }\rangle ={\displaystyle \frac{g}{\langle {\mathit{\theta }}_v\rangle }}\langle w^{\prime }{\mathit{\theta }}_v^{\prime }\rangle ,\end{array}$$

where θ_v is the virtual potential temperature and g stands for the gravity acceleration.

Turbulence production due to shear P is an important part of the budget close to the Earth's surface. Buoyancy B plays a dominant role during daytime convection; finally, the role of turbulent transport T is to transfer the kinetic energy produced close to the surface to higher altitudes.

At the beginning of the evening transition, the buoyancy flux B=0 (Nadeau et al., 2011) and the convective boundary layer collapse rapidly. Turbulence is still produced by the shear P in the surface layer; however we can assume that at higher altitudes, turbulence starts to decay freely such that the time derivative of k on the left-hand side (LHS) of Eq. (8) is balanced mainly by the dissipation ϵ. Under such assumptions, Eq. (8) rewritten in terms of the velocity scale ${\mathcal{U}}^{\mathrm{2}}=\mathrm{2}/\mathrm{3}k$ is

$$\begin{array}{}\text{(10)}& {\displaystyle \frac{\text{d}{\mathcal{U}}^{\mathrm{2}}}{\text{d}t}}=-{\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}\mathit{ϵ}\left(t\right).\end{array}$$

During the equilibrium decay, the dissipation rate is described by Taylor's law (Eq. 1) and Kolmogorov's type of turbulence kinetic energy spectrum (Eq. 3). For further comparisons with experimental data, it is convenient to express the rate of change of velocity statistics as a function of the turbulence Reynolds number, $Re=\mathcal{U}\mathcal{L}/\mathit{\nu }$. It is of note that the product of velocity and length scales 𝒰ℒ is proportional to the eddy viscosity ν_T. Hence, the Reynolds number Re in fact expresses the ratio of the eddy viscosity and the molecular viscosity $Re\sim {\mathit{\nu }}_T/\mathit{\nu }$.

After substituting Eq. (1) into the right-hand side (RHS) of Eq. (10) and further rearrangements, we obtain

$$\begin{array}{}\text{(11)}& {\displaystyle \frac{\text{d}{\mathcal{U}}^{-\mathrm{2}}}{\text{d}t}}={\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}{C}_{\mathit{ϵ}}{\displaystyle \frac{\mathrm{1}}{\mathit{\nu }}}{\displaystyle \frac{\mathit{\nu }}{\mathcal{U}\mathcal{L}}}={\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}{C}_{\mathit{ϵ}}{\displaystyle \frac{\mathrm{1}}{\mathit{\nu }}}{\displaystyle \frac{\mathrm{1}}{Re}}.\end{array}$$

To derive the corresponding equation for the rate of change in L, the equilibrium law (Eq. 1) is differentiated over time:

$$\begin{array}{}\text{(12)}& {\displaystyle \frac{\mathrm{1}}{\mathit{ϵ}}}{\displaystyle \frac{\text{d}\mathit{ϵ}}{\text{d}t}}=\mathrm{3}{\displaystyle \frac{\mathrm{1}}{\mathcal{U}}}{\displaystyle \frac{\text{d}\mathcal{U}}{\text{d}t}}-{\displaystyle \frac{\mathrm{1}}{\mathcal{L}}}{\displaystyle \frac{\text{d}\mathcal{L}}{\text{d}t}}.\end{array}$$

To express the derivative $\text{d}\mathit{ϵ}/\text{d}t$, the predictions of the classical k−ϵ turbulence model can be used (Launder and Sharma, 1974):

$$\begin{array}{}\text{(13)}& {\displaystyle \frac{\mathrm{1}}{\mathit{ϵ}}}{\displaystyle \frac{\text{d}\mathit{ϵ}}{\text{d}t}}=C_{\mathrm{0}}{\displaystyle \frac{\mathrm{1}}{{\mathcal{U}}^{\mathrm{2}}}}{\displaystyle \frac{\text{d}{\mathcal{U}}^{\mathrm{2}}}{\text{d}t}},\end{array}$$

where C₀=1.9 is a model constant. Substituting Eq. (13) into Eq. (12), we obtain

$$\begin{array}{}\text{(14)}& {\displaystyle \frac{\mathrm{1}}{\mathcal{L}}}{\displaystyle \frac{\text{d}\mathcal{L}}{\text{d}t}}=\left({\displaystyle \frac{\mathrm{3}}{\mathrm{2}}}-C_{\mathrm{0}}\right){\displaystyle \frac{\mathrm{1}}{{\mathcal{U}}^{\mathrm{2}}}}{\displaystyle \frac{\text{d}{\mathcal{U}}^{\mathrm{2}}}{\text{d}t}}.\end{array}$$

After multiplying both sides by 2ℒ², using Eqs. (1) and (10), the following relation is derived:

$$\begin{array}{}\text{(15)}& {\displaystyle \frac{\mathrm{1}}{\mathit{\nu }}}{\displaystyle \frac{\text{d}{\mathcal{L}}^{\mathrm{2}}}{\text{d}t}}=A_e{\displaystyle \frac{\mathcal{U}\mathcal{L}}{\mathit{\nu }}}=A_{e}Re,\end{array}$$

where $A_e=\mathrm{4}/\mathrm{3}(C_{\mathrm{0}}-\mathrm{3}/\mathrm{2})C_{\mathit{ϵ}}$. Hence, the classical theory predicts that the integral length scale increases with time during the decay of turbulence and that the decay of ℒ² slows down when the Reynolds number Re decreases.

Even though ℒ is expected to increase in time, the product 𝒰ℒ and also the turbulence Reynolds number Re are expected to decrease in time. Indeed, after rearranging and combining Eqs. (11) and (15), we obtain the following formula:

$$\begin{array}{}\text{(16)}& \begin{array}{rl}& {\displaystyle \frac{\text{d}\mathcal{U}\mathcal{L}}{\text{d}t}}=\mathcal{U}{\displaystyle \frac{\text{d}\mathcal{L}}{\text{d}t}}+\mathcal{L}{\displaystyle \frac{\text{d}\mathcal{U}}{\text{d}t}}=-{\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}{C}_{\mathit{ϵ}}{\mathcal{U}}^{\mathrm{2}}+{\displaystyle \frac{A_e}{\mathrm{2}}}{\mathcal{U}}^{\mathrm{2}}\\ & =-{\displaystyle \frac{\mathrm{1}}{\mathrm{15}}}{C}_{\mathit{ϵ}}{\mathcal{U}}^{\mathrm{2}}<\mathrm{0}.\end{array}\end{array}$$

2.3 Non-equilibrium decay

The scenario of non-equilibrium decay predicts that the dissipation rate scales according to the relation of Eq. (2). After introducing this relation into Eq. (10) and further rearrangements, we obtain

$$\begin{array}{}\text{(17)}& {\displaystyle \frac{\text{d}{\mathcal{U}}^{-\mathrm{2}}}{\text{d}t}}={\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}{C}_{\mathrm{ne}}{\mathcal{U}}_{\mathrm{0}}{\mathcal{L}}_{\mathrm{0}}{\displaystyle \frac{\mathrm{1}}{{\mathit{\nu }}^{\mathrm{2}}}}{\left({\displaystyle \frac{\mathit{\nu }}{\mathcal{U}\mathcal{L}}}\right)}^{\mathrm{2}}={\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}{C}_{\mathrm{ne}}R{e}_{\mathrm{0}}{\displaystyle \frac{\mathrm{1}}{\mathit{\nu }}}{\displaystyle \frac{\mathrm{1}}{R{e}^{\mathrm{2}}}},\end{array}$$

where $R{e}_{\mathrm{0}}={\mathcal{U}}_{\mathrm{0}}{\mathcal{L}}_{\mathrm{0}}/\mathit{\nu }$. Equation (17) is different than the corresponding Eq. (11) derived for the equilibrium decay. It scales with Re⁻² instead of Re⁻¹, and, additionally, it depends on the initial conditions through Re₀. During the turbulence decay, the Reynolds number Re decreases with time, hence the ratio $R{e}_{\mathrm{0}}/Re\ge \mathrm{1}$. This implies that decay rates predicted by Eq. (17) increase in time more sharply than those predicted by its equilibrium counterpart (Eq. 11). Fast anomalous changes in turbulence kinetic energy during non-equilibrium decay, which gradually decrease at later times, were observed experimentally by Meldi and Sagaut (2018).

As argued by Goto and Vassilicos (2016), in the initial stages of decay, after the forcing is stopped and long-range correlations suddenly disappear, the integral length scale starts to decrease with time. Using the self-similar form of the spectra (Eq. 5) and relation (Eq. 2), Steiros (2022b) derived the following formula for the time derivative of the length scale:

$$\begin{array}{}\text{(18)}& {\displaystyle \frac{\mathrm{1}}{\mathit{\nu }}}{\displaystyle \frac{\text{d}{\mathcal{L}}^{\mathrm{2}}}{\text{d}t}}=A_{\mathrm{ne}}-B_{\mathrm{ne}}Re,\end{array}$$

where A_ne and B_ne are coefficients related to the initial conditions and integrated spectral function. We note that the formula presented in the original paper by Steiros (2022b) was written in terms of the Reynolds number based on the Taylor length scale. However, in the non-equilibrium decay, this length scale becomes proportional to the integral length scale ℒ, which follows directly from Eq. (2); see also the discussion by Vassilicos (2015). Equation (18) is qualitatively different than its equilibrium counterpart (Eq. 15). For large turbulence Reynolds numbers Re, the time derivative of ℒ² can be negative, causing ℒ to decay in time. As the Reynolds number decreases during the decay of turbulence, the RHS of Eq. (18) will eventually become positive and ℒ will start to increase with time until $\text{d}{\mathcal{L}}^{\mathrm{2}}/\text{d}t$ reaches its maximum. At this point, the system arrives at its equilibrium state and the statistics further follow the equilibrium Taylor relation (Eq. 1), albeit with a larger value of C_ϵ≈1. This also implies that ℒ further grows in time according to Eq. (15).

Results of the numerical experiment presented in Steiros (2022b) confirmed that during the non-equilibrium decay the time derivative of ℒ² was a decreasing function of the turbulence Reynolds number. Moreover, in the initial stages of decay, $\text{d}{\mathcal{L}}^{\mathrm{2}}/\text{d}t<\mathrm{0}$.

2.4 Detection of non-equilibrium decay in the ABL

Our purpose is to show, based on the experimental evidence, that the non-equilibrium form of decay is present in the atmospheric turbulence before sunset. During non-equilibrium decay, estimating the dissipation rate from the low-wavenumber part of the wind velocity spectra (Eq. 3) or with the use of Eq. (1) becomes questionable and leads to underpredictions of the dissipation rate. On the other hand, the resolution of the Doppler lidar is not sufficient to measure small turbulent motions, which are less affected by the non-equilibrium correction. Hence, unlike in laboratory experiments, direct verification of Eqs. (1) and (2) is not possible. For this reason, the non-equilibrium will be detected indirectly, by recording changes in the scaling of the frequency spectra and the structure functions. According to Eqs. (7) and (6), increasing deviations from the Kolmogorov scaling can be explained by the decrease in the integral length scale. This is in contrast to the theory of equilibrium decay where the integral length scale should increase with time according to Eq. (15) and the scaling of spectra and structure functions should become closer to the Kolmogorov one. We calculated both frequency spectra and structure functions from time series of the vertical velocity component measured by the Doppler lidar and estimated the scaling exponents (slopes) using least-squares fitting.

Apart from non-equilibrium correction, the slopes can be affected by insufficient resolution in time and space and high noise-to-signal ratios (Frehlich, 1994; Frehlich and Cornman, 1999). Banakh et al. (2021) investigated modifications of spectra due to instrumental noise, aliasing and space averaging. In particular, in the Fourier space, the latter modification affects the whole range of scales and not only the highest wavenumbers. Moreover, the modification depends on the horizontal wind speed $\stackrel{\mathrm{‾}}{U}$. In order to convert time coordinates to space coordinates, Taylor's frozen-eddy hypothesis is used, with $x=\stackrel{\mathrm{‾}}{U}t$. This relation is justified only if the turbulence intensity defined as $\mathcal{U}/\stackrel{\mathrm{‾}}{U}$, where $\stackrel{\mathrm{‾}}{U}$ is the mean horizontal wind speed, is small enough. As $\stackrel{\mathrm{‾}}{U}$ decreases, frequency spectra become more affected.

Because of these possible modifications of the spectra, in this work we focus on detecting changes in the slopes rather than on their exact values. We filter out data with high instrumental noise and data where $\mathcal{U}/\stackrel{\mathrm{‾}}{U}>\mathrm{0.15}$. Although the frequency spectra and structure functions are mathematically equivalent, they may respond differently to different sources of errors. In parallel to the slopes, we present the calculated integral length scales and mean wind velocity to verify whether the changes in the scaling are due to changes in the integral length scale, as predicted by Eqs. (7) and (6) or, rather, are affected by changes in the mean wind speed.

3 Experimental sites, meteorological conditions and instrumentation

The data used for this work were obtained using the Doppler lidar system in rural and urban environments. The measurements were performed over a rural environment during the POLIMOS (Technical assistance for Polish Radar and Lidar Mobile Observation System) campaign, which took place between May and September 2018 at the PolWET peatland site in Rzecin (52°45^′ N, 16°18^′ E; 54 m a.s.l.), Poland. The measurements over an urban environment were performed at Warsaw Observatory Station in the center of Warsaw (52°12^′ N, 20.°58^′ E; 112 m a.s.l.), Poland. Both sites are part of the Aerosol, Clouds and Trace Gases Research Infrastructure (ACTRIS ERIC). The locations of each site are presented in Fig. 1.

During the measurements, the meteorological conditions were represented by hot and dry periods for each of the sites. In 2018, the meteorological conditions in Rzecin deviated from the reference values for precipitation and mean air temperature (512 mm and 8.63 °C, respectively), with recorded values of 464 mm and 9.63 °C (Poczta et al., 2023). Furthermore, the summer of 2018 was one of the hottest and driest periods over recent years (122 mm and 19.21 °C) in comparison to the reference values (192 mm and 18.0 °C) for this season. The meteorological conditions in 2023 in Warsaw also differed from the reference values for precipitation and mean air temperature (549.7 mm and 9.00 °C respectively), with recorded values of 620.9 mm and 11.07 °C. Even though the total amount of precipitation for 2023 in Warsaw was higher than the reference value, the summer of 2023 was drier and hotter (184.4 mm and 20.29 °C) than the reference values (257.1 mm and 17.65 °C) for this season. A relative increase in temperature (+14.96 % and +6.72 % in comparison to the reference values for Warsaw and Rzecin, respectively) and relative decrease in precipitation amount (−28.28 % and −36.46 % in comparison to the reference values for Warsaw and Rzecin, respectively) were observed for both of the measurement locations during the summer seasons. It follows that the meteorological conditions for both stations were relatively similar and were characterized by higher temperatures and less precipitation for the summer season for both 2018 and 2023 in comparison to the historical data. The reference values for precipitation and mean air temperature for Rzecin were calculated based on the Szamotuły-Baborówko station meteorological data (IMGW-PIB, 2024a) for the period between 1990–2014 (further data not available). The reference values for Warsaw were calculated from “Climate Standards 1990–2020” (IMGW-PIB, 2024b). The amount of precipitation and average air temperature value for the summer season of 2023 in Warsaw were calculated based on the Warsaw-Filtry station meteorological data (IMGW-PIB, 2024a).

Figure 1 Location of ACTRIS ERIC experimental sites, the Rzecin PolWET station of the Poznań University of Life Sciences and Warsaw Observatory Station of the University of Warsaw. The maps were provided by the National Geoportal (https://www.geoportal.gov.pl/, last access: 20 March 2024). The authors acknowledge Zuzanna Rykowska (University of Warsaw) for preparing the figure.

For each location, the vertical and horizontal wind profiles were obtained using measurements from StreamLine (HALO Photonics) Doppler lidars. In the rural environment, the measurements were taken using the Doppler lidar provided by the Atmospheric Physics Group of the University of Granada (GFAT-UGR). The Doppler lidar operating in Warsaw is owned by the Remote Sensing Laboratory (RS-Lab) at the Faculty of Physics of the University of Warsaw. The lidars comprise a solid-state pulsed laser emitting at 1.5 µm and a heterodyne detector with fiber-optic technology. The emission is provided with pulses of energy at 100 µJ; a pulse duration of 200 ns; and pulse repetition rates of 15 and 10 kHz for the GFAT-UGR and RS-Lab systems, respectively. The signal acquisition is performed in continuous and autonomous vertical mode, and regular measurements are made in the vertical azimuth display (VAD). A more detailed description of the Doppler lidar system can be found in Pearson et al. (2009).

The lidar signal acquisition was performed continuously in vertical mode, obtaining a vertical wind component with 30 m spatial and 1 s temporal resolution. For the horizontal profiles of the wind, the VAD scans with a constant elevation of 70° and 12 azimuth points were performed every 30 min. The focus of the optical system was at 535±35 m (Ortiz-Amezcua et al., 2022) and at the infinite value for GFAT-UGR and RS-Lab lidars, respectively. To support our findings, we additionally analyzed the momentum and heat fluxes (30 min averages). The eddy covariance observations were performed in parallel to the lidar measurements with instruments mounted on a meteorological tower (4.5 $\mathrm{m}\mathrm{a}.\mathrm{g}.\mathrm{l}.$) at the Rzecin PolWET station (see Poczta et al., 2023, for details) and on the Radiation Transfer Laboratory measuring platform on the roof of the building of the Faculty of Physics in Warsaw (22 $\mathrm{m}\mathrm{a}.\mathrm{g}.\mathrm{l}.$).

4 Methodology

The whole database of Doppler lidar and surface flux measurements consists of 4 months (June–September 2018) of measurements in a rural environment and 4 months (June–September 2023) of measurements in an urban environment. To obtain the data ready for further analyses, the raw Doppler lidar data were firstly background-corrected, using the calibration procedures proposed by Manninen et al. (2016) and Vakkari et al. (2019), and secondly filtered out by values of the signal-to-noise ratio threshold (SNR = 1.006) (Manninen et al., 2016). Errors in the vertical and horizontal velocity measurements were calculated using the software-processing HALO lidar toolbox (Manninen, 2019), with methods proposed by Rye and Hardesty (1993) and Pearson et al. (2009) for the vertical components and by Päschke et al. (2015) for the horizontal ones.

The lidar data were first used to estimate the ABL height. We compare the results of two methods. The first method is the gradient method of the backscatter signal. In lidar measurements, the backscatter signal in the ABL is significantly stronger than in the free troposphere, so there is a distinct change in backscatter signal values when it passes through the boundary between the ABL and the free troposphere (Wang et al., 2020; Liu et al., 2022). The boundary layer height is defined as the height where the range-corrected signal sharply decreases, that is, where its gradient has a minimum value (Dang et al., 2019). The second variance method was used to estimate the height of the convective boundary layer (Dewani et al., 2023). The top of the layer is estimated as the height where the variance of the vertical velocity decreases below a certain threshold value. We use the threshold of 0.09, as in Dewani et al. (2023).

To compute the slopes of the frequency spectra and structure functions and values of integral length scales, the vertical velocity measurements at different heights were grouped in 0.5 h intervals. The recorded signal was decomposed into the mean and fluctuating part as

$$\begin{array}{}\text{(19)}& w^{\prime }=w-\langle w\rangle ,\end{array}$$

where 〈w〉 is a 600 s running average. This detrending removes the largest convective scales from outside the inertial sub-range and allows for a better convergence of statistics within 0.5 h intervals. The turbulence velocity scale was calculated from the vertical velocity fluctuations as

$$\begin{array}{}\text{(20)}& \mathcal{U}=\langle {w^{\prime }}^{\mathrm{2}}{\rangle }^{\mathrm{1}/\mathrm{2}}\end{array}$$

for each height and time interval. Cases where the turbulence intensity was larger than 0.15, which would not support Taylor's hypothesis, were filtered out. The frequency spectra were computed for each group, and a logarithmic fit in the frequency range $f\in [\mathrm{0.15}{\mathrm{s}}^{-\mathrm{1}}$, 0.3 s⁻¹] was performed without fixing the slope. In the logarithmic plot, the power-law function forms a straight line, and its slope is equal to the scaling exponent. We calculated the slopes at each height and each time interval using the least-squares algorithm and investigated whether the slopes deviate from Kolmogorov's predictions ($-\mathrm{5}/\mathrm{3}$ and $\mathrm{2}/\mathrm{3}$ scaling of the frequency spectrum and structure function, respectively). To determine the integral length scale, we first calculated the two-point transverse correlation coefficient:

$$\begin{array}{}\text{(21)}& g\left(r\right)={\displaystyle \frac{\langle w\left(x\right)w(x+r)\rangle }{\langle w(x{)}^{\mathrm{2}}\rangle }}.\end{array}$$

According to the theory of homogeneous, isotropic turbulence (Pope, 2000), the function g(r) takes the form

$$\begin{array}{}\text{(22)}& g\left(r\right)=\exp \left(-{\displaystyle \frac{r}{\mathcal{L}}}\right)\left(\mathrm{1}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{r}{\mathcal{L}}}\right)\end{array}$$

and crosses the horizontal axis at r=2ℒ. We numerically integrated the function g(r) from r=0 to r=2ℒ. This integral should be approximately equal to 0.57ℒ.

Figure 2 presents the methodology to obtain the mean horizontal wind, slopes of the frequency spectrum and structure function, and the integral length scale for each interval. For slopes, we additionally use the R² threshold (O'Connor et al., 2010) such that fits with R² less than 0.6 are considered noise and discarded.

To calculate the buoyancy flux $\langle w^{\prime }{b}^{\prime }\rangle$ defined in Eq. (9) using the data from the PolWET station, we approximated the virtual potential temperature as

$$\begin{array}{}\text{(23)}& {\displaystyle }{\mathit{\theta }}_{\mathrm{v}}=\mathit{\theta }\left[\mathrm{1}+\left({\displaystyle \frac{R_{\mathrm{v}}}{R_{\mathrm{d}}}}-\mathrm{1}\right)q_{\mathrm{v}}\right]\approx \mathit{\theta }\left[\mathrm{1}+\mathrm{0.61}{q}_{\mathrm{v}}\right],\end{array}$$

where R_v and R_d are gas constants for water vapor and dry air, respectively, and q_v is the mixing ratio of water vapor; we assumed θ≈T, where T is absolute temperature. With this the buoyancy flux was calculated as

$$\begin{array}{}\text{(24)}& \langle w^{\prime }{b}^{\prime }\rangle ={\displaystyle \frac{g}{\langle {\mathit{\theta }}_v\rangle }}\langle w^{\prime }{T}^{\prime }\rangle +\mathrm{0.61}g\langle w^{\prime }{q}_{\mathrm{v}}^{\prime }\rangle .\end{array}$$

For the Warsaw data, we approximated the virtual potential temperature with the measured sonic temperature θ_v≈T_s and calculated the buoyancy flux as

$$\begin{array}{}\text{(25)}& \langle w^{\prime }{b}^{\prime }\rangle ={\displaystyle \frac{g}{\langle T_{\mathrm{s}}\rangle }}\langle w^{\prime }{T}_{\mathrm{s}}^{\prime }\rangle .\end{array}$$

Figure 2 Methodology chain of obtaining turbulence properties: calibration (dark yellow), filtering procedures (orange), calculating turbulence properties (blue) and remaining procedures (pale yellow).

5 Results

5.1 ABL height, velocity and timescales

As we focus on the evening hours in this work, we denote the sunset time as t=0. Turbulence statistics are calculated as 0.5 h averages or medians relative to the sunset time for both Rzecin and Warsaw for the same months during the summer. We present the statistics as functions of time to analyze how they change during rapid decay of turbulence.

We first present the evolution of the heights of the ABL and convective boundary layer (CBL) in Fig. 3. The uncertainties were estimated as a standard error of the mean (SEM) of the half-hour time intervals; see Appendix A. As is observed, the ABL height at the rural site is lower and decreases more rapidly compared to the urban site. This may indicate the influence of the urban heat island effect in large cities. The behavior of the CBL height is similar; i.e., much higher values are observed in Warsaw. In both the Rzecin and the Warsaw sites, the CBL collapses rapidly ca. 2.5 h before sunset. The results are in agreement with the long-term study by Wang et al. (2020).

Figure 3 The median values of the ABL and CBL heights and error estimates (see Appendix A) during the summer season (June–September) at the Rzecin PolWET station in 2018 (a) and Warsaw Observatory Station in 2023 (b).

In order to estimate the beginning of the evening transition, we present the buoyancy fluxes, their median values and standard errors in Fig. 4. At $t=-\mathrm{2}$ h, the median values of the fluxes are still positive. Afterwards they cross the zero level at the Rzecin site and become close to zero at the Warsaw site.

Figure 4 The median values of the buoyancy flux and error estimates (see Appendix A) during the summer season (June–September) at the Rzecin PolWET station in 2018 (a) and Warsaw Observatory Station in 2023 (b).

The buoyancy fluxes and the CBL height are further used to calculate the characteristic scales which govern the flow during daytime convection. The convective Deardorff scale w_* and the corresponding timescale τ are defined as (Deardorff, 1970)

$$\begin{array}{}\text{(26)}& w_{*}={\left(D\langle w^{\prime }{b}^{\prime }{\rangle }_{\mathrm{0}}\right)}^{\mathrm{1}/\mathrm{3}},\mathit{\tau }={\displaystyle \frac{D}{w_{*}}},\end{array}$$

where the CBL height is denoted by D and $\langle w^{\prime }{b}^{\prime }{\rangle }_{\mathrm{0}}$ is the surface value of the buoyancy flux. We present both scales in Fig. 5 and additionally compare w_* with the friction velocity u_*. A description of the error estimates is given in Appendix A. As is seen in Fig. 5, both in Rzecin and in Warsaw, w_* is still larger than u_* at $t=-\mathrm{2}$ h and the timescale τ increases sharply around $t=-\mathrm{2}$ h. This timescale is larger in Warsaw, which suggests that the turbulence decay is slower in the urban surface layer. Shortly after $t=-\mathrm{2}$ h, u_* becomes larger than or comparable to w_* (within the confidence intervals), which implies that turbulence production due to shear P becomes dominant in the surface layer. We will assume that at higher altitudes the contribution of the shear production is negligible and that, to the leading order, the evolution of the turbulence kinetic energy is described by Eq. (10). Under such assumptions, the changes in 𝒰 and ℒ should be described by formulas derived in Sect. 2.2 and 2.3. For this reason we will further treat statistics at $t=-\mathrm{2}$ h as initial conditions and focus our attention on the time interval of ±2 h relative to the sunset.

Figure 5 The mean friction velocity u_*, convective Deardorff scale w_*, timescale τ and error estimates (see Appendix A) during the summer season (June–September) at the Rzecin PolWET station in 2018 (a, b) and Warsaw Observatory Station in 2023 (c, d).

5.2 Spectra and structure functions

At 18:30 UTC, with convection still present, the wind energy spectrum was measured at a height of 195 $\mathrm{m}\mathrm{a}.\mathrm{g}.\mathrm{l}.$, i.e., within the mixed layer (see Fig. 6a). In this case the spectrum is visibly steeper than $-\mathrm{5}/\mathrm{3}$. The steep slopes of frequency spectra in the convective regime were also reported by other authors (Darbieu et al., 2015). Xie and Huang (2022) speculated that these steep slopes could be linked to the presence of inverse cascades at large scales, which leads to the $-\mathrm{11}/\mathrm{5}$ Bolgiano–Obukhov scaling (Bolgiano, 1959). However, as argued by Banakh et al. (2021), in the case of poorly resolved data, the steepening of slopes may also be caused by the artificial instrumental dissipation due to effective low-pass filtering.

Figure 6b shows that before sunset, when the convective layer rapidly decays, the frequency spectrum becomes less steep than the Kolmogorov's prediction. This observation is not related to instrumental artificial dissipation, which rather causes the opposite effect (steepening of the spectra).

Figure 6 Exemplary frequency spectra of vertical wind measured on 28 June 2018 at the Rzecin PolWET station, at (a) 195 $\mathrm{m}\mathrm{a}.\mathrm{g}.\mathrm{l}.$ at 18:30 UTC and (b) 195 $\mathrm{m}\mathrm{a}.\mathrm{g}.\mathrm{l}.$ at 19:30 UTC, showing the steeper and less steep frequency spectrum, respectively, as compared to Kolmogorov's prediction.

In Fig. 7, the time–height evolution of the slopes of the frequency spectra and velocity structure functions and the estimate of the CBL height are shown. The regions marked by red and orange colors in Fig. 7a and yellow and light green in Fig. 7b exhibit scaling steeper than that of Kolmogorov. Before sunset, when the CBL collapses, the turbulence kinetic energy decays rapidly (El Guernaoui et al., 2019) and a sharp decrease in the slopes at all heights is observed. Decrease in the slopes is also seen in the upper part of the convective ABL, where the rising updrafts become weaker. Therein, stable stratification possibly alters the spectra and structure functions. The stratification effects on spectral slopes were included in a recent model by Cheng et al. (2020). This issue is, however, beyond the scope of the present work, as we focus instead on the modification of the spectra due to nonstationarity.

Figure 7 Slopes of the frequency spectra and velocity structure functions and estimates of the CBL height for data measured on 28 June 2018 at the Rzecin PolWET station. The boxes marked in black indicate the slopes of frequency spectra equal to $-\mathrm{5}/\mathrm{3}\pm \mathrm{1}/\mathrm{6}$ and slopes of the velocity structure function equal to $\mathrm{2}/\mathrm{3}\pm \mathrm{1}/\mathrm{12}$.

5.3 Influence of the detrending window

We next investigated how different values of the detrending window influence the results in order to choose the most appropriate value. The size of the detrending window can potentially affect the turbulence velocity scale and the integral length scale, calculated from Eqs. (20) and (22), respectively. To estimate the errors, we took into account the standard error of the mean and errors in velocity estimates from the HALO lidar toolbox (see Appendix A for details). We calculated and compared both mean and median values, and their comparison is shown in Fig. 8. Median values are advantageous if the probability distribution of the data is non-Gaussian or in the presence of rare but very large or very small values (outliers), which affect the mean. As can be seen in Fig. 8, all the mean values are larger then the corresponding medians. Both the mean and the median values of 𝒰 are not much affected by the change in the detrending window in the algorithm. They decrease with time during the evening transition. Mean and median values of the integral length scales increase with increasing size of the detrending window. The data become considerably scattered for the largest window. This is to be expected, as ℒ converges slowly and a large sample is needed to estimate it with sufficient accuracy (Lenschow et al., 1994). The time span needed to calculate statistics is proportional to ℒ; hence larger length scales require longer time spans. As a result, 0.5 h averages may be insufficient to reduce random errors to acceptable levels needed to calculate the time derivatives from Eqs. (15) and (18). Independent of the size of the detrending window, we find that the median values of integral length scales first decrease with time and next increase or become constant. This conclusion seems to be universal and allows further analysis to be performed using only one set value for the detrending window, which has been chosen to be 600 s. The median values of the product 𝒰ℒ (and also the turbulence Reynolds number Re) decrease with time, as expected, for all detrending windows. The corresponding mean values decrease in time before sunset for the two smaller detrending windows.

Figure 8 Mean and median values of the turbulence velocity scale 𝒰, the integral length scale ℒ and the product 𝒰ℒ with error estimates (see Appendix A) for different detrending windows for data measured at the Rzecin PolWET station in August 2018.

5.4 Turbulence properties: statistical analysis

To investigate how the turbulence properties change not only in time but also with altitude, we divided results into three altitude ranges: 75–345 $\mathrm{m}\mathrm{a}.\mathrm{g}.\mathrm{l}.$, 345–585 $\mathrm{m}\mathrm{a}.\mathrm{g}.\mathrm{l}.$ and 585–855 $\mathrm{m}\mathrm{a}.\mathrm{g}.\mathrm{l}.$ According to Fig. 3, the measurements are within the residual layer, and in Rzecin the level 75–345 $\mathrm{m}\mathrm{a}.\mathrm{g}.\mathrm{l}.$ corresponds to the lower and middle part of the ABL, 345–585 $\mathrm{m}.\mathrm{a}.\mathrm{g}.\mathrm{l}.$ reaches the top of the layer, and the third level of 585–855 $\mathrm{m}.\mathrm{a}.\mathrm{g}.\mathrm{l}.$ is placed partly above the mean top of the ABL. The ABL is higher in Warsaw such that the three levels corresponds to its lower, middle and upper part, respectively. In spite of the scatter of results for individual 0.5 h averages, median values of the slopes of frequency spectra determined for the range $f\in [\mathrm{0.15}{\mathrm{s}}^{-\mathrm{1}}$, 0.3 s⁻¹] (Fig. 9) clearly increase from values close to Kolmogorov's $-\mathrm{5}/\mathrm{3}$ to around −1 after sunset. Analogously, slopes of the structure functions (Fig. 10) determined for the range r∈[40 m, 150 m] decrease with time for both sites, Rzecin and Warsaw. The inertial ranges of structure functions tend to be smaller than those for wind frequency spectra. Hence, due to the finite temporal resolution of the measurements, the calculated structure functions may be affected by large eddies from beyond the inertial range. As a result, the slopes of structure functions are gentler than the Kolmogorov $\mathrm{2}/\mathrm{3}$ predictions even at $t=-\mathrm{2}\mathrm{h}$, especially for the Rzecin site. The slopes also show the altitude-dependent relation. What can be observed from both figures (Figs. 9 and 10) is that for Warsaw the median values of slopes at higher altitudes are much closer to Kolmogorov's predictions than they are for the rural case. This can be partly explained by the urban heat island effect and differences between the urban and rural morphologies. The turbulent flows are generated due to an intense shear at the top of the canopy layer H (Roth, 2000). In the urban environment, this is at the height of the rooftops of the buildings. In Warsaw, in a 1 km × 1 km area centered on the RS-Lab, H≈12.2 m. In Rzecin, the top canopy layer height reaches H=2–2.5 m. As far as the thermally driven turbulent flows are concerned, the city centers are usually warmer than the suburbs and areas of urban agglomeration because of the urban heat island phenomenon (Kuchcik et al., 2014; Stopa-Boryczka et al., 2002) and, related to this, have relatively high heat capacity (Oke, 1982) compared to the rural environment. This implies that the urban surface can still emit heat even after sunset. Heat emitted from a warm urban surface generates convection and mixes the air in the urban canopy layer. It also generates a dome of warm air in higher parts of the boundary layer. The temperature profile of this dome is quasi-adiabatic and similar to the temperature profiles around midday (Oke, 1987). Results presented in Figs. 9 and 10 suggest, however, that the heat island effects are most significant 2 h before sunset at the highest altitude range of 585–855 $\mathrm{m}\mathrm{a}.\mathrm{g}.\mathrm{l}.$ Afterwards, a rapid change in the scaling exponent is observed.

Figure 9 Slopes of the frequency spectra of vertical wind, with means, medians and error estimates (see Appendix A) at different altitudes for data measured during the summer season (June–September) in rural and urban environments at the Rzecin PolWET station in 2018 (a, c, e) and Warsaw Observatory Station in 2023 (b, d, f), respectively.

Figure 10 Slopes of the second-order structure functions, means, medians and error estimates (see Appendix A) at different altitudes for data measured during the summer season (June–September) in rural and urban environments at the Rzecin PolWET station in 2018 (a, c, e) and Warsaw Observatory Station in 2023 (b, d, f), respectively.

In order to examine if the changes in the scaling exponent are induced by changes in the mean velocity, we analyzed how $\stackrel{\mathrm{‾}}{U}$ changes with time and height. If turbulence were Kolmogorov-like and the spectra were affected only by spurious modifications due to insufficient resolutions in time and space, then a decrease in mean velocity would increase the absolute value of the slopes by introducing artificial dissipation. An increase in the mean velocity, on the other hand, would bring the scaling closer to the Kolmogorov $-\mathrm{5}/\mathrm{3}$ (or $\mathrm{2}/\mathrm{3}$ for the structure function) but not below this value. Figure 11 presents the mean horizontal velocity $\stackrel{\mathrm{‾}}{U}$. As is seen in the figure, $\stackrel{\mathrm{‾}}{U}$ increases with time mostly close to the surface and only after sunset. At higher altitudes, it has only a slight tendency to increase. At the same time, the absolute values of slopes in Figs. 9 and 10 decrease considerably below $-\mathrm{5}/\mathrm{3}$ and $\mathrm{2}/\mathrm{3}$, respectively. Hence, we conclude that changes in the slopes are not primarily affected by the changes in the mean wind speed. Instead, the collapse of the largest convective motions possibly leads to the non-equilibrium states of turbulence as predicted by Eqs. (6) and (7). Later on, the size of the inertial range decreases and is shifted towards small scales (large wavenumber) which are not detected by the lidars.

A difference between rural and urban sites observed in Fig. 11 is that in the latter, mean velocity starts to slightly increase before sunset at low altitudes. As described by Mahrt (2017), when turbulence decreases rapidly, the airflow becomes more influenced by the surface heterogeneity and horizontal temperature variations. The temperature variations lead to local horizontal pressure gradients. These, in turn, could induce stronger horizontal winds to increase the horizontal temperature transport.

Mahrt (1981) discussed the formation of early-evening calm periods during which the mean velocity decreased in the surface layer and had a tendency to increase at higher altitudes. Busse and Knupp (2012), on the other hand, observed a decrease in the mean speed at altitudes of up to 500 m. The early-evening calm periods in the surface layer were also recorded in other studies (Mahrt et al., 2012; Román-Cascón et al., 2015). No systematic decrease in the mean wind speed is observed in Fig. 11 before sunset. However, we recall that, to calculate turbulence properties, we removed data for which the Taylor frozen-eddy hypothesis was not satisfied, that is data with a high $\mathcal{U}/\stackrel{\mathrm{‾}}{U}$ ratio. This procedure could, to some extent, affect the tendencies observed in Fig. 11.

Figure 11 Horizontal velocity, with means, medians and error estimates (see Appendix A) at different altitudes for data measured during the summer season (June–September) in rural and urban environments at the Rzecin PolWET station in 2018 (a, c, e) and Warsaw Observatory Station in 2023 (b, d, f), respectively.

The dependence of the integral length scales on time and altitude is presented in Fig. 12. The vertical integral length scales increase with the altitude. Moreover, at the urban site they are larger at higher altitudes, which follows from stronger convection and more shear-driven turbulent flows in this part of the ABL. As argued by Akinlabi et al. (2022), the roughness surface layer in cities may be higher than previously expected and can reach up to $z/H=\mathrm{30}$, where H is the mean height of buildings (at the Warsaw site, H≈12.2 m).

Figure 12 Integral length scales at different altitudes for data, with means, medians and error estimates (see Appendix A), measured during the summer season (June–September) in rural and urban environments at the Rzecin PolWET station in 2018 (a, c, e) and Warsaw Observatory Station in 2023 (b, d, f), respectively.

According to Fig. 13, the turbulent velocity scale decreases with time for both rural and urban cases. In connection with the stable horizontal wind velocity, before and after sunset the figure shows that the observed turbulent flows are mostly driven by the convection, which agrees with our observational statement about the decay of the ABL and the presence of the decaying convection-driven turbulent flows before sunset.

Figure 13 Standard deviation of the vertical velocity component at different altitudes for data, with means, medians and error estimates (see Appendix A), measured during the summer season (June–September) in rural and urban environments at the Rzecin PolWET station in 2018 (a, c, e) and Warsaw Observatory Station in 2023 (b, d, f), respectively.

To study decay of turbulence with reference to the values measured 2 h before sunset, normalized 15 min median values of $\mathcal{L}/{\mathcal{L}}_{\mathrm{0}}$ and $\mathcal{U}/{\mathcal{U}}_{\mathrm{0}}$, where ${\mathcal{L}}_{\mathrm{0}}=\mathcal{L}(t=-\mathrm{2}\mathrm{h})$ and ${\mathcal{U}}_{\mathrm{0}}=\mathcal{U}(t=-\mathrm{2}$ h), are plotted in Fig. 14. In this figure we include only data measured before sunset. It is seen that both turbulence kinetic energy and the integral length scale have a tendency to decrease with time before sunset. As follows from Eqs. (6) and (7), a decrease in the integral length scale during the decay of turbulence will cause increasing deviations from the Kolmogorov scaling. This appears to be the case in the observed decay of the convective boundary layer and explains the changes in the slopes, presented in Figs. 9 and 8.

Figure 14 Integral length scales and turbulence velocity scales before sunset normalized by the values measured at $t=-\mathrm{2}\mathrm{h}$, with means, medians and error estimates (see Appendix A): the Rzecin PolWET station, rural environment (a, c, e), and Warsaw Observatory Station, urban environment (b, d, f).

We finally calculated the time derivatives $\text{d}{\mathcal{U}}^{-\mathrm{2}}/\text{d}t$ and $\text{d}{\mathcal{L}}^{\mathrm{2}}/\text{d}t$ from the 0.5 h median values (for both sites and at all heights) for −2 h $\le t\le \mathrm{0}$ (evening hours) and present them in Fig. 15 as functions of Re. We compare experimental data with the classical predictions of Eqs. (11) and (15) with C_ϵ=0.5 and the non-equilibrium predictions (17) and (18). The coefficients A_ne, B_ne and C_ne in Eqs. (17) and (18) were estimated from linear regression as the best fit.

The errors in the estimates in Fig. 15 are considerable such that it is difficult to identify the scaling of $\text{d}{\mathcal{U}}^{-\mathrm{2}}/\text{d}t$. However, $\text{d}{\mathcal{U}}^{-\mathrm{2}}/\text{d}t$ is clearly a decreasing function of the turbulence Reynolds number Re, as predicted by both equilibrium and non-equilibrium relations.

As far as $\text{d}{\mathcal{L}}^{\mathrm{2}}/\text{d}t$ is concerned, points follow the non-equilibrium scaling. As the Reynolds number decreases, both equilibrium and non-equilibrium predictions for $\text{d}{\mathcal{L}}^{\mathrm{2}}/\text{d}t$ become close, so we expect that at even lower Re values, the data will follow the equilibrium predictions. Hence, by analysis of both panels in Fig. 15, it can be concluded that the non-equilibrium decay is likely to be present in the initial stages.

Figure 15 (a) Dependence of $\text{d}{\mathcal{U}}^{-\mathrm{2}}/\text{d}t$ on the Reynolds number. (b) Dependence of $\text{d}{\mathcal{L}}^{\mathrm{2}}/\text{d}t$ on the Reynolds number. Derivatives are estimated from half-hour medians and error estimates (see Appendix A) (symbols with error bars) compared to equilibrium relations (Eqs. 11 and 15) (solid lines) and non-equilibrium relations (Eqs. 17 and 18) (dotted lines).

5.5 Dissipation rates: comparison

Due to large deviations from the equilibrium scaling, the energy dissipation estimates with the use of Eq. (3) or Eq. (4) are not reliable at $t>-\mathrm{2}$ h. Therefore we used Eq. (4) to calculate the profile of ϵ only at $t=-\mathrm{2}$ h and for $z/D<\mathrm{0.6}$. An example of the second-order structure function is presented in Fig. 16a, and the dissipation rate, non-dimensionalized with the use of w_* and the CBL height D, is presented in Fig. 16b and c, as a function of $z/D$. We use these estimates for comparison with the predictions of Eqs. (1) and (2) in Fig. 17 and mark them as red dots at $t=-\mathrm{2}$ h. Figure 17 presents median dissipation values before sunset obtained using the classical Taylor law (Eq. 1), with C_ϵ=0.5 and assuming the non-equilibrium scenario. In the latter, the dissipation rate is expected to follow the non-equilibrium relation (Eq. 2), until its values become equal to the predictions of the classical Taylor law (Eq. 1) but with a higher value of the constant C_ϵ=1, which is the upper bound of the dissipation coefficient (Bos and Rubinstein, 2018). The dissipation rate further follows the equilibrium Taylor law (Eq. 1) with C_ϵ=1. During this latter period, the frequency spectra will still deviate from Kolmogorov's scaling at low wavenumbers (Steiros, 2022b). The same is true for the structure functions; see Eq. (7). The equilibrium $-\mathrm{5}/\mathrm{3}$ and $\mathrm{2}/\mathrm{3}$ slopes of spectra and structure functions will be reached asymptotically at high wavenumbers which are not measured by the lidar system. By comparing both scaling laws, it can be concluded that the equilibrium Taylor law underpredicts the dissipation rate of turbulence kinetic energy by up to a factor of 2.

The dissipation values decrease with altitude and in time. In Rzecin, the decrease is most rapid closer to the surface. Under the conditions of radiative cooling over a rural area, surface-based inversion is developed, which stops the air mixing above the ground faster than at the urban sites. As a consequence of very similar meteorological conditions for both stations (dry and warm conditions), we assume that the differences in turbulence properties result mainly from differences in surface morphology between both environments. Due to the difference in surface morphology, at the urban site the effect of wind shear and the presence of shear-driven turbulent flows are still significant at altitudes of 75–345 $\mathrm{m}\mathrm{a}.\mathrm{g}.\mathrm{l}.$ On the other hand, Fig. 17 suggests that at higher altitudes, the dissipation rates decrease in time much faster over the urban site. Shortly before sunset, the heat island effects are possibly much more limited in height.

Figure 16 (a) Exemplary second-order structure function at $t=-\mathrm{2}\mathrm{h}$ and (b, c) the non-dimensionalized turbulence kinetic energy dissipation rates ${\mathit{ϵ}}^{+}=\mathit{ϵ}D/w_{*}^{\mathrm{3}}$ at $t=-\mathrm{2}\mathrm{h}$ estimated from Eq. (4) (solid line) with uncertainty ranges (dashed lines) – see Appendix A.

Figure 17 The turbulence kinetic energy dissipation rate calculated from the equilibrium Taylor law (Eq. 1) and with the non-equilibrium scenario, with error estimates (see Appendix A) at different altitudes for data measured in rural and urban environments at the Rzecin PolWET station (a, c, e) and Warsaw Observatory Station (b, d, f), respectively. Dots at $t=-\mathrm{2}$ h are estimates from the second-order structure function.

6 Conclusions

The analysis performed in this work showed time, altitude and surface-type dependencies of the properties of turbulent flows during the evening transition at a time span of ±2 h relative to the sunset. The calculated slopes of the frequency spectra and structure function of vertical wind before sunset deviate increasingly from Kolmogorov's predictions, which agrees with the non-equilibrium scenario. These values also deviate more with altitude, implying the possible maximum height of the turbulence presence in the ABL during its rapid decay. We argued that it is possible to explain the increasing deviations of the slopes with the use of recent theories of turbulence. The crucial part was the observed decrease in the turbulence length scale, ℒ, during turbulence decay, which was also predicted by the non-equilibrium relations.

In this work we assumed that within the residual layer, where turbulence decays, the production and turbulent transport of kinetic energy are negligible such that, to the leading order, the time change in kinetic energy is balanced by the dissipation rate. Under these assumptions, the rates of change $\text{d}{\mathcal{U}}^{-\mathrm{2}}/\text{d}t$ and $\text{d}{\mathcal{L}}^{-\mathrm{2}}/\text{d}t$ can be expressed as functions of the turbulence Reynolds number Re. We choose as the initial condition $t=-\mathrm{2}$ h, when the estimated convective Deardorff scale w_* is still larger than the friction velocity u_*. For this choice, our results suggest that statistics follow non-equilibrium relations before sunset, when the turbulence Reynolds numbers are very high. Hence, non-equilibrium relations should be taken into account for the estimation of the dissipation rate of turbulence kinetic energy in the initial stages of decay.

Our work also shows differences in the turbulence properties between two different environments, the rural and urban. The latter is much more morphologically diverse and has a higher heat capacity. We found that over the urban area, turbulence is initially present at higher altitudes. Moreover, the convective timescale τ calculated 2 h before sunset is larger in Warsaw, which suggests that the decay of turbulence in the urban layer is slower as compared to the rural one. However, our observations suggest that this is the case only at lower altitudes. At heights of 585–855 $\mathrm{m}\mathrm{a}.\mathrm{g}.\mathrm{l}.$, turbulence seems to decrease very rapidly over the urban site. Turbulence production by shear, which is affected by the surface morphology, could also contribute to the difference between the two environments, especially at lower altitudes. We conclude that the influence of the surface morphology and heat capacity on turbulence decay is significant and should be accounted for in parametrization schemes, which is in line with the results reported by Couvreux et al. (2016).

Appendix A: Error estimates

In the analyses we used the standard error of the mean (SEM) of a quantity 〈ϕ〉, defined as

$$\begin{array}{}\text{(A1)}& \mathrm{\Delta }\langle \mathit{\varphi }{\rangle }_{\mathrm{sem}}={\displaystyle \frac{\mathit{\sigma }}{\sqrt{N}}},\end{array}$$

where σ is the standard deviation and N is the number of samples. We assumed that the relative errors in median values are equal to the relative errors of the mean values. The uncertainties presented in Figs. 4–17 were determined as follows:

• Errors in the ABL and CBL heights D presented in Fig. 3 were estimated as SEMs.

• Errors in the buoyancy fluxes from Fig. 4 were determined from a sum of relative SEMs of temperature and temperature fluxes as

  $$\begin{array}{}\text{(A2)}& {\displaystyle \frac{\mathrm{\Delta }\langle w^{\prime }{b}^{\prime }\rangle }{\langle w^{\prime }{b}^{\prime }\rangle }}={\displaystyle \frac{\mathrm{\Delta }\langle {\mathit{\theta }}_v{\rangle }_{\mathrm{sem}}}{\langle {\mathit{\theta }}_v\rangle }}+{\displaystyle \frac{\mathrm{\Delta }\langle w^{\prime }{\mathit{\theta }}_v^{\prime }{\rangle }_{\mathrm{sem}}}{\langle w^{\prime }{\mathit{\theta }}_v^{\prime }\rangle }}.\end{array}$$

• Error in the friction velocity u_* in Fig. 5 was estimated as SEMs; errors in the Deardorff velocity w_* and the timescale τ are defined as

  $$\begin{array}{}\text{(A3)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{\Delta }{w}_{*}}{w_{*}}}={\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}\left({\displaystyle \frac{\mathrm{\Delta }\langle w^{\prime }{b}^{\prime }\rangle }{\langle w^{\prime }{b}^{\prime }\rangle }}+{\displaystyle \frac{\mathrm{\Delta }{D}_{\mathrm{sem}}}{D}}\right),\\ & {\displaystyle \frac{\mathrm{\Delta }\mathit{\tau }}{\mathit{\tau }}}={\displaystyle \frac{\mathrm{\Delta }{w}_{*}}{w_{*}}}+{\displaystyle \frac{\mathrm{\Delta }{D}_{\mathrm{sem}}}{D}}.\end{array}\end{array}$$

• To estimate errors in the turbulence velocity and length scales in Figs. 8, 12, 13 and 14, we took into account errors in vertical velocity measurements δw obtained from the HALO lidar toolbox (Manninen, 2019) and errors of the mean. The half-hour averages ${\mathcal{U}}_{\mathrm{1}/\mathrm{2}}^{\mathrm{2}}$ were calculated by averaging over time, and their measurement error $\mathit{\delta }{\mathcal{U}}_{\mathrm{1}/\mathrm{2}}^{\mathrm{2}}$ is determined as

  $$\begin{array}{rl}& {\mathcal{U}}_{\mathrm{1}/\mathrm{2}}^{\mathrm{2}}\pm \mathit{\delta }{\mathcal{U}}_{\mathrm{1}/\mathrm{2}}^{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }t}}\underset{t}{\overset{t+\mathrm{\Delta }t}{\int }}{\left[w\left(t\right)+\mathit{\delta }w\left(t\right)\right]}^{\mathrm{2}}\text{d}t\\ & \approx {\mathcal{U}}_{\mathrm{1}/\mathrm{2}}^{\mathrm{2}}\pm {\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }t}}\underset{t}{\overset{t+\mathrm{\Delta }t}{\int }}\mathrm{2}w\left(t\right)\mathit{\delta }w\left(t\right)\text{d}t,\end{array}$$

  where Δt=0.5 h. Next, the mean values were calculated by averaging over N half-hour averages. Error of the mean value was calculated from the sum of squares of the measurement error δ𝒰² and the convergence error $\mathrm{\Delta }{\mathcal{U}}_{\mathrm{sem}}^{\mathrm{2}}$ as

  $$\begin{array}{}\text{(A4)}& \begin{array}{rl}& \mathrm{\Delta }{\mathcal{U}}^{\mathrm{2}}={\left[{\left(\mathit{\delta }{\mathcal{U}}^{\mathrm{2}}\right)}^{\mathrm{2}}+{\left(\mathrm{\Delta }{\mathcal{U}}_{\mathrm{sem}}^{\mathrm{2}}\right)}^{\mathrm{2}}\right]}^{\mathrm{1}/\mathrm{2}},\\ & \mathit{\delta }{\mathcal{U}}^{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^N\mathit{\delta }{\mathcal{U}}_{\mathrm{1}/\mathrm{2}}^{\left(i\right)\mathrm{2}}.\text{Further,}\\ & {\displaystyle \frac{\mathrm{\Delta }\mathcal{U}}{\mathcal{U}}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{\mathrm{\Delta }{\mathcal{U}}^{\mathrm{2}}}{{\mathcal{U}}^{\mathrm{2}}}}.\end{array}\end{array}$$

  The length scales ℒ were also calculated based on the vertical velocity measurements, and we assumed that their error consists of the measurement error δℒ and the convergence error Δℒ_sem as

  $$\begin{array}{}\text{(A5)}& \begin{array}{rl}& \mathrm{\Delta }\mathcal{L}={\left[{\left(\mathit{\delta }\mathcal{L}\right)}^{\mathrm{2}}+{\left(\mathrm{\Delta }{\mathcal{L}}_{\mathrm{sem}}\right)}^{\mathrm{2}}\right]}^{\mathrm{1}/\mathrm{2}},\\ & \text{where}{\displaystyle \frac{\mathit{\delta }\mathcal{L}}{\mathcal{L}}}={\displaystyle \frac{\mathit{\delta }{\mathcal{U}}^{\mathrm{2}}}{{\mathcal{U}}^{\mathrm{2}}}}.\end{array}\end{array}$$

• Errors in the slopes in Figs. 9 and 10 and errors in the horizontal velocity in Fig. 11 were calculated as SEMs.

• Derivatives presented in Fig. 15 were estimated as forward finite differences, and their errors were estimated as

  $$\begin{array}{}\text{(A6)}& \begin{array}{rl}& \mathrm{\Delta }{\displaystyle \frac{\text{d}{\mathcal{L}}^{\mathrm{2}}}{\text{d}t}}=\mathrm{2}{\displaystyle \frac{\mathcal{L}(t+\mathrm{\Delta }t)\mathrm{\Delta }\mathcal{L}(t+\mathrm{\Delta }t)}{\mathrm{\Delta }t}}\\ & +\mathrm{2}{\displaystyle \frac{\mathcal{L}\left(t\right)\mathrm{\Delta }\mathcal{L}\left(t\right)}{\mathrm{\Delta }t}},\\ & \mathrm{\Delta }{\displaystyle \frac{\text{d}{\mathcal{U}}^{-\mathrm{2}}}{\text{d}t}}=\mathrm{2}{\displaystyle \frac{{\mathcal{U}}^{-\mathrm{3}}(t+\mathrm{\Delta }t)\mathrm{\Delta }\mathcal{U}(t+\mathrm{\Delta }t)}{\mathrm{\Delta }t}}\\ & +\mathrm{2}{\displaystyle \frac{{\mathcal{U}}^{-\mathrm{3}}\left(t\right)\mathrm{\Delta }\mathcal{U}\left(t\right)}{\mathrm{\Delta }t}}.\end{array}\end{array}$$

• Errors in the dissipation rate estimates in Fig. 16 are SEMs.

• Errors in the dissipation rates presented in Fig. 17 are

  $$\begin{array}{}\text{(A7)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{\Delta }\mathit{ϵ}}{\mathit{ϵ}}}=\mathrm{3}{\displaystyle \frac{\mathrm{\Delta }\mathcal{U}}{\mathcal{U}}}+{\displaystyle \frac{\mathrm{\Delta }\mathcal{L}}{\mathcal{L}}}\text{and}\\ & {\displaystyle \frac{\mathrm{\Delta }\mathit{ϵ}}{\mathit{ϵ}}}=\mathrm{2}{\displaystyle \frac{\mathrm{\Delta }\mathcal{U}}{\mathcal{U}}}+\mathrm{2}{\displaystyle \frac{\mathrm{\Delta }\mathcal{L}}{\mathcal{L}}}\end{array}\end{array}$$

  in the equilibrium and non-equilibrium cases, respectively.