Asymptotic limits of the attached eddy model derived from an adiabatic atmosphere

Table 1 summarizes the formulations of $\sigma_{u}$ and $\sigma_{u}^{2}$ found from field experiments in the ASL. According to Monin-Obkhov Similarity Theory (MOST), the streamwise velocity standard deviations can be expressed as [43, 39]

where $L$ is the Obukhov length measuring the height at which mechanical production of turbulent kinetic energy (TKE) balances buoyancy production or destruction of TKE, and $\phi_{u}$ is the so-called stability correction function that varies with $(z-d)/L$. Under neutral conditions where $|L|\rightarrow\infty$, $\phi_{u}$ reduces to

Several ASL measurements conducted under near-neutral conditions as well as laboratory flows have been used to fit the empirical parameter $C_{0}$. [33] surveyed several experiments and stated that $C_{0}$ appears independent of $z$ but varies with terrain and with values ranging between 2.1 and 2.45. The survey also noted that under varying stability conditions, vertical and horizontal velocity components exhibit distinct behaviours. While the vertical velocity standard deviation aligns well with MOST predictions, streamwise velocity components are affected by $z-d$ and by $1/L$ differently. Specifically, under near-convective conditions, a change in $z$ has a negligible effect on streamwise velocity standard deviations, whereas a change in $L$ has a pronounced effect. A recent study attributed this limitation of MOST to the anisotropy of the Reynolds stress tensor [54].

To account for the ‘non-MOST’ behaviour of streamwise velocity components and the increase of streamwise wind fluctuations with the deepening of boundary layer height under near-convective conditions, [44] proposed an empirical formulation based on ASL observations from 30 m to 87 m given by

This formulation became widely cited in the atmospheric science literature [62] although it assumed that changes in $z$ have a negligible effect on $\sigma_{u}$. It is to be noted that under neutral conditions Eq. 6 is similar to Eq. 5. Another extensive study was undertaken for unstable conditions [61] where TKE production and TKE dissipation rates were in approximate balance. These experiments reported a mild $z$-dependence of the streamwise velocity variance and a refinement to Panofsky’s formulation was proposed given as

Those developments were viewed as adjustments to MOST and made no apparent contact with the AEM. In fact, after the work of [23], [8] commented that MOST formulations for the ASL appear to have missed predictions from the AEM about the role of $\delta$ in $\sigma_{u}^{2}/u_{*}^{2}$ [7], which is briefly reviewed next.

The AEM postulates that turbulence in wall-bounded flows can be described by a hierarchy of self-similar eddies attached to the wall [56]. These eddies are geometrically similar at different scales with their size proportional to their distance from the wall. Using the equilibrium-layer hypothesis where the mechanical production of TKE balances the viscous dissipation, the AEM predicts a specific scaling law in the energy spectrum [47, 48]. Further, by integrating the streamwise velocity spectrum across all wavenumbers, $\sigma_{u}^{2}$ can be formulated as in Eq. 2 for the overlap region of turbulent flows. This relation is a cornerstone of the AEM, linking the spectral characteristics of turbulence at a given $z$ to the profile of $\sigma_{u}^{2}$ [56, 46].

Since the classic book by [56], the AEM has drawn significant experimental and theoretical interest (see Table 1). [47] expanded the AEM to include near-wall regions, addressing inner flow dynamics and providing a theoretical and experimental foundation for the logarithmic law. Subsequently, [49] introduced a $C{Re^{+}}^{-1/2}$ term to Eq 2 to account for different surface types—smooth, rough, and wavy—identifying distinct coefficients ($A_{1}$, $B_{1}$, and $C$) for each. In a further advancement, [48] incorporated a viscous correction term $V[Re^{+}]$, arguing for the formulation’s independence from Reynolds number variations. Building on these foundations, [36] proposed a similarity relation to describe the streamwise velocity variance in the logarithmic region, considering inviscid attached eddies and incorporating a viscous correction term $V_{g}[Re^{+},\frac{z}{\delta}]$ in the inner region and a wake correction term $W_{g}[\frac{z}{\delta}]$ in the outer region. This formulation was evaluated using wind tunnel and near-neutral ASL data with success.

[37] tested this similarity formulation using data collected at the Surface Layer Turbulence and Environmental Science Test (SLTEST) facility, confirming the logarithmic behaviour of $\sigma_{u}^{2}$ in near-neutral ASL conditions. [41] further explored this formulation by neglecting the correction terms for the viscous and outer flow effects, which were deemed insignificant in their study. These experiments matched well with AEM predictions when setting $A_{1}=1.03$, consistent with the prior value reported in [48].

Further experimental support for Eq. 2 was provided by [35], who considered four datasets including boundary layers, pipe flow and ASL measurements with $2\times 10^{4}<Re_{\tau}<6\times 10^{5}$. Their results not only affirmed the presence of a universal logarithmic region but also estimated $A_{1}$ to be around 1.26 for lab flow and 1.33 for ASL data, which is different from the previous estimation of $A_{1}=1.03$. [52] reported an $A_{1}=1.26$ based on wind tunnel data within the range $6000<Re_{\tau}<20000$. Similarly, [58] found that their ASL experiments supported the estimate of $A_{1}=1.33$, while [19] used SLTEST data near the wall (just above the roughness layer) to study the high-order moments of the streamwise velocity, finding a fitted $A_{1}$ of about 0.9. This value appears to be closer to the value reported by [25] for the near-neutral ASL derived from streamwise velocity spectra.

Despite the variability in $A_{1}$, the debate continues regarding the origin of the log-law of the streamwise velocity variance. Some studies have employed spectral budget analysis and dimensional analysis, complemented by laboratory experiments [42, 4], to explain the logarithmic behaviour of streamwise velocity variance. Recently, [21] introduced a model inspired by Townsend’s original work, addressing the spectrum in the logarithmic layer for various $z/\delta$ values. This model suggests that the coefficients $A_{1}$ and $B_{1}$ depend on the Reynolds number. Their analysis of wind tunnel data at $Re_{\tau}\sim 6123\text{–}19680$ indicated that the approximated $A_{1}$ values vary between 1.01 and 1.09. Furthermore, their model challenges the notion that a ‘-1’ power law in the streamwise energy spectrum is necessary for the logarithmic behaviour in streamwise velocity variance. This evolving understanding underscores the need to verify the universality of the logarithmic behaviour of streamwise velocity variance, particularly within the ASL and assess to what extent the coefficient $A_{1}$ can be derived from the spectrum of the streamwise velocity at a single level $z$, considered next.

Assuming the sonic anemometers are nearly level, requiring only minor tilt corrections to align the vertical axis perpendicular to the mean flow, the double-rotation method [60] involves two sequential rotations: first, a pitch rotation to set the mean vertical velocity component ($\overline{w}$) to zero, and second, a yaw rotation to set the mean lateral velocity component ($\overline{v}$) to zero. Given the relatively flat and homogeneous terrain, this standard method was initially applied to correct the tilt of the anemometers.

Since achieving perfectly levelled sonic anemometers on the 62-meter tower is challenging, the planar fit method [60] was applied to assess the consistency of turbulent quantities between the two coordinate rotation methods. This method, which is more robust to initial misalignments, fits a plane to the wind data over an extended period, based on specific wind direction sectors.

In the analysis of ASL data, fluctuations with periods less than about one hour are generally considered turbulence, while slower fluctuations are synoptic-scale and treated as part of the mean flow [63]. Consequently, a one-hour averaging period was selected as a compromise between the need to resolve large eddies reliably and stationarity considerations. However, [37] suggests that the neutral periods are often short, on the order of several minutes. Previous studies have used various averaging periods depending on the stability and conditions of the boundary layer. For instance, 10 minutes were used in the stable boundary layer in the CASES-99 experiment [6], 30 minutes in the unstable conditions when mechanical production of TKE was compared to viscous dissipation [61], 20 minutes in neutral conditions in the SLTEST experiment [37], and 1 hour in neutral conditions in the same SLTEST experiment [20]. Additionally, 30 minutes were used in SLTEST [19], 15 minutes in the unstable LATEX experiment [32], and 1-minute for stable and 30 minutes for unstable conditions across 13 datasets [54].

A shorter sampling period, while increasing steadiness in mean meteorological conditions, can also affect statistical convergence by reducing sample size, increasing sensitivity to outliers, and potentially biasing the representative nature of the data. In general, the difference between ensemble-averaging (the sought quantity) and time-averaging is labelled as systematic bias. This bias declines with reduced $2T_{L}/T_{p}$ [31], where $T_{L}$ is the integral time scale of a flow variable and $T_{p}$ is the averaging period. Typical $T_{L}$ values for the streamwise velocity are on the order of 1-2 minutes and hence selecting $T_{p}=60$ minutes is acceptable for reducing the systematic bias. However, potential contamination by trends and changes in the meteorological conditions to the near-neutral ASL turbulence is acknowledged, and some trend removal must still be performed to remove synoptic weather phenomena.

After the coordinate rotation, the time series in Figure 2a shows a superposition of the high-frequency turbulent fluctuations and low-frequency oscillations. De-trending is employed to remove long-term trends or patterns in the data, which can result from various factors such as instrument drift and atmospheric changes [38]. By removing the low-frequency content, the analysis can focus on short-term variability (i.e., turbulence) and ensure a certain amount of stationarity.

A common method in micrometeorology is linear detrending, where the line of best fit over the 1-hour averaging period (dashed white line in Figure 2b) is subtracted from the original time series. While primarily affecting the low-frequency part of the signal, linear de-trending can impact all frequencies and introduce oscillations at higher frequencies [38].

Alternatively, high-pass filtering can isolate low-frequency signals by convolving the original data with a transfer function in the frequency domain, a method often used in post-processing [40, 20, 37]. However, determining a clear cut-off frequency to separate turbulent motions from long-term trends is challenging, posing a risk of filtering out low-wavenumber turbulent motions. A $10^{th}$ order Butterworth filter was applied in this work with a cut-off frequency equivalent to a 2000-m cut-off wavelength. It translates to a filter size of 100 to 350 seconds that varies with mean wind speed at different levels. Those sizes align closely with the observed time period of low-frequency signals in the time series (indicated by the two vertical lines in Figure 2b). This filter size is consistent with previous turbulent boundary layer studies that selected 180 seconds in [20] and 181.8 seconds in [51].

As shown in Figure 2c, the linear detrended streamwise velocity energy spectrum $S_{u}(f)$ levels off as $f\to 0$. In contrast (see Figure 2d), the high-pass filtered $S_{u}(f)$ shows more removal of low-frequency content (very large scale motions) than linear de-trending, with a much faster decay as $f\to 0$ as expected. Nevertheless, the best method depends on site conditions. Therefore, both methods are initially applied, and a comparison is performed to ensure the derived AEM coefficients are not sensitive to the de-trending approach.

After de-trending, hourly averaging of all instantaneous data are conducted to obtain first the mean component and then the turbulent fluctuations as $u^{\prime}=u-\overline{u}$, $v^{\prime}=v-\overline{v}$, $w^{\prime}=w-\overline{w}$, and $T^{\prime}=T-\overline{T}$, where $T$ is the air temperature. The Reynolds stress methods and profile methods are two conventional approaches to estimating the friction velocity in the ASL and those two methods will also be compared for reference.

1. 1.

   Reynolds stress method In an ideal ASL that is high Reynolds number, stationary, planar homogeneous, lacking subsidence, and characterized by a zero mean pressure gradient, the momentum fluxes must be invariant to variations in $z$. To test the constant flux layer assumption and determine $u_{*}$ to be used as a normalizing velocity in the AEM, a local friction velocity is defined as

   $$u_{*}=[\overline{u^{\prime}w^{\prime}}^{2}+\overline{v^{\prime}w^{\prime}}^{2}]^{1/4}.$$

   (11)

2. 2.

   Profile method Another approach to estimating $u_{*}$ is based on the logarithmic law of the mean velocity as shown in Eq. 1. Measurements of the mean velocity profile can be used to estimate the friction velocity by fitting the measured $\overline{u}$ against $\ln(z)$. To do so, it is necessary to identify the beginning and ending points that follow a log-linear relation, which can involve some subjectivity. The magnitude of the friction velocity is directly related to the choice of the von Kármán constant [30].

Agreement between these two methods is also used when identifying the span of the ASL.

With the $u_{*}$ and $\sigma_{u}^{2}$ within the ISL measured, a linear regression is conducted to derive the values of $A_{1}$ by fitting the $\sigma_{u}^{2}$ to $\ln(z)$ utilizing Eq. 2. Notably, this fitting method for obtaining $A_{1}$ does not require direct measurement of $\delta$ but cannot be used to infer the $B_{1}$ constant of the AEM. Figure 6a demonstrates one such linear regression for the hour of 1500LT on Sep. 25, 2020, detrended by high pass filtering. The logarithmic mean velocity profile is also shown for comparison. The results indicate that the fitted $\kappa$ (=0.41) is very close to the accepted 0.4 value. Moreover, the $\sigma_{u}^{2}$ profile aligns with AEM predictions confirming the logarithmic dependence on $z$ of streamwise velocity variance within the ISL (highlighted by closed circles) at a very high Reynolds number.

To estimate $B_{2}$ in Eq. 3, $\sigma_{w}^{2}$ is fitted against the local $u_{*}^{2}$ within the ISL, and the resulting slope is retrieved. If the vertical velocity energy spectrum is approximated by two regions: a flat region from $k=0$ to $k=k_{a}$ (energy splashing) and an inertial subrange region from $k=k_{a}$ to $k=\infty$, then the area under this idealized spectrum with $k_{a}=(\kappa z)^{-1}$ and TKE production balancing TKE dissipation is $(\sigma_{w}/u_{*})^{2}=(5/2)C_{o,w}$. Here, $C_{o,w}$ is the Kolmogorov constant for the one-dimensional vertical velocity spectrum in the inertial subrange (=0.65). This $B_{2}$ equals 1.65 (instead of 1.39 inferred here, see Figure 6b) is thus anticipated for infinite Reynolds number when the transition wavenumber $k_{a}=(\kappa z)^{-1}$ and the energy splashing region due to wall effects extends all the way from $k_{a}$ to $k=0$. In the area calculations of the vertical velocity spectrum, it is interesting to note that $(3/2)C_{o,w}$ originate from extending the inertial subrange from $k=k_{a}$ to $k=\infty$ whereas only $C_{o,w}$ is the contribution from the large eddies (i.e. eddies exceeding $\kappa z$). Thus, a plausibility argument to the reduced $B_{2}$ measured here is that the 10 Hz sampling frequency and anemometer path averaging may disproportionately under-resolve the inertial subrange components of $w$ compared to its longitudinal velocity counterpart. Nonetheless, similar values for $B_{2}$ (i.e., around 1.4) have been observed in ASL literature [64].

To assess the consistency of the logarithmic behaviour of the streamwise velocity and associated AEM coefficients in the near-neutral ASL, the same fitting procedures are applied across all 120 neutral cases. The comparison of the results between linearly detrended data and high-pass filtered data can be found in Table 3. The analysis shows that the estimated values of $\kappa$ for both methods are similar, yet smaller than 0.4, aligning with previous experiments [1, 34] as well as other theories based on the co-spectral budget [29] where $\kappa$ was linked to the Rotta constant, an isotropization constant linked to the rapid distortion theory, and the Kolmogorov constant. Furthermore, the fitted $B_{2}$ are relatively consistent across two detrending methods. As a result, the focus of the current analysis is on $A_{1}$. The high-pass filtered data exhibit a more robust logarithmic scaling of $\sigma_{u}^{2}$ than the linearly detrended data. This outcome is based on the finding that the linear detrending method yields a wider range of $R^{2}$ values with an average of 0.67 and a standard deviation of 0.32. In contrast, the high-pass filtering shows both a higher mean (0.81) and a reduced standard deviation (0.23) in $R^{2}$ scores. Moreover, the derived $A_{1}$ values of linear-detrended data exhibit significant variability, often departing considerably from the conventional range of 1 to 1.3, and frequently resulting in negative $A_{1}$ values (31.63% and 9.18% negative occurrence for linear detrending and high-pass filtering, respectively). The potential reasons for this deviation are discussed next.

Given an averaging period of 1 hour, the turbulence data might experience unsteadiness associated with larger-scale meteorological influences such as static pressure gradients. According to [18], ‘inactive’ turbulence becomes more significant as the pressure gradient increases, potentially obscuring the scaling behaviour predicted by the AEM. While the linear-detrending method only removes the mean trend, it might not sufficiently eliminate the (non-linear) large-scale meteorological influences, potentially leading to negative fitted values of $A_{1}$. This underscores the importance of removing larger-scale meteorological influences in ASL data when comparing them to a canonical turbulent boundary layer. Furthermore, sensitivity tests indicate that the fitting outcomes are not sensitive to the choice of the cut-off frequency in the 100-200 second range for the high-pass filtering method. Consequently, the following analysis focuses on high-pass filtered data.

Finally, Table 3 presents results from both the double rotation and planar fit methods, showing no substantial differences between them except for the fitted $A_{1}$, which may be more sensitive to coordinate transformations due to its higher moment nature. As shall be seen in Section III.2.2, once data quality control is applied, the difference between the double rotation and planar fit methods becomes much small even for fitted $A_{1}$. For simplicity, the following discussion is primarily based on the double rotation method.

To understand the variations in the fitted $A_{1}$, key factors that influence data quality are investigated including statistical fitting performance, non-stationarity effects, the presence of a constant $\sigma_{w}$ with height in the ISL, and the Reynolds number effect.

1. 1.

   Fitting Performance Before discussing the fitted coefficients of the AEM to the data here, it is necessary to evaluate the goodness of the logarithmic scaling for $\sigma_{u}^{2}$. The coefficient of determination ($R^{2}$) from the logarithmic fitting is examined and anomalous cases with $R^{2}$ less than 0.6 are filtered out to ensure that only high-quality fits are considered when evaluating $A_{1}$.

2. 2.

   Non-Stationarity Effects Stationarity is a necessary condition for the AEM and for the identification of the ISL. Two indices of stationarity (IST) are defined here as:

   $$\mathrm{IST_{wspd}}=\frac{|\frac{1}{12}\sum_{i=1}^{12}{\sigma_{u}^{i}}^{2}-{\sigma_{u}}^{2}|}{{\sigma_{u}}^{2}}$$

   (12)

   $$\mathrm{IST_{wdir}(\degree)}=\mathrm{|wdir_{5min}-wdir_{1hr}|}_{max}$$

   (13)

   In Eq. 12, ${\sigma_{u}^{i}}^{2}$ is the streamwise velocity variance for each 5-min block within the 1-hour run and ${\sigma_{u}}^{2}$ is the variance for the entire hour. In Eq. 13, the notation $|...|_{max}$ represents the maximum difference between the wind direction observed at a 5-minute interval and the wind direction observed at a 1-hour interval among the 12 individual blocks. When the value of $\mathrm{IST}$ is small (close to zero), the streamwise velocity time series can be viewed as stationary. As the value of $\mathrm{IST}$ increases, the non-stationarity effects become significant. According to [38], $\mathrm{IST_{wspd}}<30\%$ should be achieved to ensure stationarity of the $u$ time series and this metric is employed here.

3. 3.

   Presence of a Constant $\sigma_{w}$ An index $I_{\sigma_{w}}=|\sigma_{w}-\left<\sigma_{w}\right>|/\left<\sigma_{w}\right>$ is defined to measure the deviation of the vertical velocity standard deviation at each level from its depth-averaged value ($\left<\sigma_{w}\right>$) across the ISL. A small value suggests that the vertical velocity standard deviation is vertically uniform consistent with the AEM. The necessary conditions for the vertical velocity variance to be invariant with $z$ in the ISL may be derived from the mean vertical velocity equation for stationary, planar homogeneous flow in the absence of subsidence at very high Reynolds number. For these idealized conditions, the vertical velocity variance is given by

   $\displaystyle\frac{\partial\sigma_{w}^{2}}{\partial z}=-\left(\frac{1}{\rho}\right)\left(\frac{\partial\overline{P}}{\partial z}\right)-g,$

   (14)

   where $\overline{P}$ is the mean pressure. When $\overline{P}=-\rho gz$ (i.e. hydrostatic), $\partial\sigma_{w}^{2}/\partial z=0$ or $\sigma_{w}^{2}$ is constant with respect to $z$ within the ISL. That is, the AEM requires $\overline{P}$ to be hydrostatic above and beyond the zero-pressure gradient condition needed to ensure a constant turbulent stress with $z$. The effects of adverse pressure gradients on $C_{1}$ have already been reported in wind tunnels and pipes [57]. In the absence of mean pressure gradients, values of $C_{1}$ varying from 0.90 to 0.92 have been reported. Interestingly, for large adverse pressure gradients, a -1 power law was still reported but values of $C_{1}$ as high as 17 were computed [57]. These laboratory experiments underscore connections between the aforementioned constant $\sigma_{w}^{2}$ with $z$, finite mean pressure gradients, and the numerical values of $C_{1}$ and $A_{1}$.

4. 4.

   Reynolds Number As earlier noted, sufficiently large Reynolds numbers are necessary when applying the AEM allowing for self-similarity of turbulent structures [56]. Numerous studies highlight the importance of high Reynolds numbers in validating the AEM and observing the expected turbulent behaviour [34, 53, 4, 19]. Therefore, the dependency of fitted $A_{1}$ on $Re_{\tau}$ is also investigated. To calculate a bulk Reynolds number analogous to what is reported in laboratory studies, $\delta$ must be estimated. No direct measurement of $\delta$ was available during the experiment and the boundary layer height was estimated from $\delta=C_{a}u_{*}/f_{c}$, where $C_{a}$ is a dimensionless empirical constant typically varying between 0.1 and 0.3. Here $C_{a}=0.1$ is used after examining the streamwise energy spectrum as discussed later. The $f_{c}$ is the Coriolis parameter and is given by $f_{c}=2\Omega\sin\left(\phi\right)$, where $\Omega$ is the angular velocity of the Earth and $\phi$ is the latitude of the location [65].

Before examining how fitted $A_{1}$ is affected by these individual factors mentioned above, it is noted that these factors are not entirely independent of each other. Figure 7a reveals a moderate negative correlation between $R^{2}$ and $\mathrm{IST_{wspd}}$, indicating that increased non-stationarity in wind speed reduces the quality of the logarithmic fit. A weak negative correlation between $R^{2}$ and $\mathrm{IST_{wdir}}$ suggests a discernible decrease in fit quality with higher wind direction non-stationarity. More dynamically interesting is that higher computed Reynolds numbers ($Re_{\tau}$) are positively correlated with a better fitting quality ($R^{2}$), agreeing with the basic assumption of AEM. Non-stationarity in wind speed ($\mathrm{IST_{wspd}}$) is negatively correlated with $Re_{\tau}$, and greater deviations in the vertical velocity standard deviation ($I_{\sigma_{w}}$) and higher non-stationarity indices tend to occur at lower Reynolds numbers. Overall, a smaller $Re_{\tau}$ reflects, to some extent, the non-stationarity effect and the absence of a constant $\sigma_{w}$, all of which implicitly affect the fitting quality for the constants of the AEM at very high $Re_{\tau}$. Therefore, the analysis primarily focuses on the role of $Re_{\tau}$ in describing the variations in fitted $A_{1}$ values.

Figure 7b illustrates that the derived $A_{1}$ values vary appreciably at lower Reynolds numbers but tend to asymptotically converge to a relatively constant range around 1.25 as ${Re}_{\tau}$ increases. This suggests that despite the ASL’s Reynolds numbers being on the order of $10^{7}$, the fitted $A_{1}$ values remain dependent on the Reynolds number. It is important to note that this dependence is not associated with scale separation due to the energy cascade as suggested by previous studies [34, 53, 4, 19], but rather results from the combined effects of non-stationarity, the presence or absence of the overlap region, and systematic bias.

When $Re_{\tau}>4\times 10^{7}$ is used as the criterion, 41 cases remain, all of which satisfy $\mathrm{IST_{wspd}}<30\%$. Among these 41 cases, 39 cases have $R^{2}>0.6$, with the 2 cases of $R^{2}<0.6$ shown as open circles. Therefore, the subsequent analysis focuses on the 39 cases where $R^{2}>0.6$ and $Re_{\tau}>4\times 10^{7}$.

Table 3 and  4 present the results before and after applying the criteria of $Re_{\tau}>4\times 10^{7}$ and $R^{2}>0.6$. After these quality controls are applied, the mean value of $A_{1}$ changes from 1.25 to 1.22. The mean value of $\kappa$ remains consistent (from 0.37 to 0.38), though with a reduced standard deviation. The $R^{2}$ values show marked improvement, and the non-stationarity indices decrease substantially. The vertical variability of $\sigma_{w}$ also decreases from 0.03 to 0.02.

Post-screening, the patterns of near-neutral occurrence hours and wind directions remain consistent with the previous data. Instances that are infrequent pre-screening become even rarer after the screening, while common instances become more frequent (see Figure 3c and 3d). The double rotation and planar fit methods show even smaller differences (see Table 3). The analysis indicates that the variations in $A_{1}$ are primarily driven by non-stationarity considerations and Reynolds number dependence. However, the variability of $A_{1}$ (std equal 0.33) remains non-trivial compared to the classical values derived from laboratory experiments despite all the data quality controls imposed on these ASL measurements. This highlights the complexity of ASL flows, necessitating case-by-case investigations, which will be discussed in Section III.3.2.

In unstable conditions, thermal plumes enhance vertical mixing and modify the turbulent eddy structure. In stable conditions, near-surface turbulence and the residual layer above it progressively decouple with increasing stability, reducing the influence of ABL-scale motions on near-surface turbulence and diminishing the size of streamwise streaks. In those conditions, near-surface turbulence becomes intermittent [12]. These modifications can influence the streamwise velocity variance and energy spectra, including the $k^{-1}$ scaling [5]. Although $\left|z/L\right|<0.1$ is applied to select near-neutral cases, the selected cases could still be influenced by buoyancy effects (especially the higher $z$ measurements). Discerning these effects is now considered.

It is observed that sensible heat flux ($|H|$) shows a strong positive correlation with $Re_{\tau}$ (correlation coefficient of 0.65), indicating that higher sensible heat flux tends to be associated with higher Reynolds numbers. After the data quality control, the sensible heat flux remains high and variable both before and after screening, posing a challenge in achieving both high Reynolds numbers and minimal buoyancy effects. In the ABL, this correlation is not surprising as stronger mean winds usually occur during daytime conditions when surface heating is large. Figure 8a illustrates the relation between the stability parameter and the fitted $A_{1}$. The screening criteria of $Re_{\tau}>4\times 10^{7}$ and $R^{2}>0.6$ exclude data with stability parameters greater than 0.026 (positive, stable) and less than -0.05 (negative, unstable). The remaining data are roughly evenly distributed across the range (-0.05, 0.025) with no significant trend. This suggests that there is little stability effect on $A_{1}$ after screening when stability is quantified using the local Obukhov length and the local measurement height. This finding is not sensitive to the initial choice of $\left|z/L\right|<0.1$, which justifies its use to indicate near-neutral conditions.

In contrast, the influence of stability conditions on the fitted values of $\kappa$ is visible in Figure 8b. Even after applying the stability constraint, the fitted $\kappa$ values remain influenced by variations in atmospheric stability despite the near-neutral filter, consistent with other prior experiments [1]. It is possible to correct the stability effect on $\kappa$ by employing MOST and the Businger-Dyer relation for the stability correction function [9]. After correcting the stability effects, the median value of $\kappa$ estimated from the dataset analyzed here is 0.36 (not shown). This aligns with Figure 8 where the vertical dashed line intersects the regression line roughly at $\kappa=0.36$.

The pre-multiplied streamwise velocity energy spectra are computed within the ISL using the Welch method [59] while applying a Hann window with a $2^{14}$ window size and $50\%$ overlap between consecutive windows. The power spectral density is then interpolated onto 500 logarithmically spaced frequencies, spanning from $10^{-4}$ to 5 Hz (Nyquist frequency). Figure 9 shows the ensemble average (i.e., average over 39 post-screening cases) of the pre-multiplied $u$-spectra. The streamwise wavenumber $k$ is calculated based on Taylor’s hypothesis as $k=2\pi f/\overline{u}$, where $f$ is the frequency. Due to the normalization of wavenumbers by $z$, the pre-multiplied spectra at different heights collapse well in the inertial subrange [24]. The pre-multiplied $u$-spectrum rises at the tail due to high-frequency noise. However, this high-frequency noise does not significantly affect the streamwise velocity variance and is neglected.

A $k^{-2/3}$ power-law scaling is evident in the inertial subrange wavenumber region, consistent with Kolmogrov’s theory. More interestingly, a plateau of about 1 to 1.25 begins roughly at $k\delta=1$ (see Figure 9a), consistent with the onset of the $k^{-1}$ scaling predicted by the AEM, suggesting that the estimation of $\delta=C_{a}u_{*}/f_{c}$ with $C_{a}=0.1$ is plausible. Furthermore, this plateau occurs within the range $\mathcal{O}(0.1)<kz<\mathcal{O}(1)$ as shown in Figure 9b, indicating a robust $k^{-1}$ power-law scaling for $E_{uu}(k)$. However, the $k^{-1}$ scaling region is just under a decade in extent prompting the question of whether individual cases show clear $k^{-1}$ scaling and whether the plateau values ($C_{1}$) correspond to the fitted $A_{1}$ values, which will be discussed next.

To analyze the linkage between the logarithmic behaviour of streamwise velocity variance and the $k^{-1}$ scaling in the streamwise energy spectra ($E_{uu}$), each of the 39 remaining cases is examined individually. $C_{1}$ are approximated by the 95^th percentile values of $kE_{uu}/u_{*}^{2}$ for simplicity. Although this is an approximation, it enables us to explore a general correlation between $A_{1}$ and $C_{1}$.

The profile of the vertical velocity variance and the TKE budget are also examined. In the ISL, the TKE budget simplifies to a balance between the mechanical production ($P$) and the viscous dissipation rate ($\epsilon$). The production rate is determined by:

where ${\partial U}/{\partial z}$ is approximated in two ways: first, using ${u_{*}}/{\kappa z}$ based on the log law, where $\kappa$ is the fitted value; second, using a third-order polynomial fitting of the mean wind profile over all 12 levels. These two methods yield similar results (as shown in Figure A1a in the appendix). The production rate obtained using the second method is used in the analysis of the TKE budget. The dissipation rate ($\epsilon$) is determined by both the second-order structure functions ($D_{2}(r)$) and the third-order structure functions ($D_{3}(r)$), respectively. At very high Reynolds numbers, these are given as follows

where $C_{2}$ is a coefficient ($\approx 1.97$ for one-dimensional wavenumber), and $r$ represents the separation distance in the streamwise direction, set to be the minimum between 1 meter and $z/2$ [10]. Prior studies have shown that the ratio of $D_{2}(r)$ in the vertical to the streamwise direction most closely matches the isotropic value of 4/3 when $0.63$ m $\leq r\leq$ 1.0 m [10]. A similar eddy size range was also identified in other ASL experiments [26]. Figure A1b shows that the dissipation rates derived from the second-order structure function are slightly smaller than their third-order counterpart. In the following analysis, both dissipation rates are used to assess the sensitivity of the AEM to TKE dissipation rate estimation.

Three cases (out of 39) exhibit anomalous $\sigma_{u}^{2}$ profiles (Group 1; see Figure A2), while nine cases show anomalous $\sigma_{w}^{2}$ profiles with clear trends within the ISL (Group 2; see Figure A3). Seven cases display anomalous $E_{uu}$ profiles, where the pre-multiplied streamwise velocity spectrum presents two or more peaks in the production range (Group 3; see Figure A4), suggesting the contribution of eddies other than the attached eddies. Five cases have an anomalous TKE balance (Group 4; see Figure A5), with the ratio of production to dissipation rates showing clear trends in the ISL, likely indicating non-equilibrium turbulence, potentially influenced by sudden changes in wind speed, temperature gradients, or other factors. These anomalous behaviours of turbulent properties can coincide, though this overlap is not fully explored in the categorization. These anomalies may be due to local topographic features, instrumentation issues, or transient meteorological phenomena. It is acknowledged that the categorization is subjective to the authors’ choices. Therefore, this section serves only to provide a qualitative conclusion.

Ultimately, nine cases are flagged as high-quality data (Group 5), serving as benchmarks for understanding the connection between streamwise velocity variance and the streamwise velocity energy spectrum. Figure 10 provides an example of this high-performance group on Mar. 29, 2021, at 0400LT. The profile of $\sigma_{u}^{2}$ exhibits a logarithmic behaviour, with a fitted $A_{1}$=1.12 (Figure 10a). Similar to the ensemble average, the pre-multiplied $E_{uu}$ reach a plateau in the range of $\mathcal{O}(0.02)<kz<\mathcal{O}(1)$. The estimated $C_{1}$ is 1.23, marginally higher than the fitted $A_{1}$. At the inertial subrange, the streamwise energy spectra align well with a $k^{-2/3}$ scaling (Figure 10b) as expected. The vertical velocity variance remains constant with $z$ in the ISL (Figure 10c), with the $\sigma_{w}^{2}/u_{*}^{2}$ values centred around 1.3. The ratio $P/\epsilon$ derived from both $D_{2}(r)$ and $D_{3}(r)$ is invariant with height within the ISL. However, $P/\epsilon$ based on $D_{2}(r)$ is smaller than that based on $D_{3}(r)$, with the latter being near-unity.

To investigate whether a $C_{1}=A_{1}$ relation can be observed in the adiabatic ASL, the estimated $A_{1}$ and $C_{1}$ of the 39 post-filtered cases were grouped based on these additional quality controls and plotted against $Re_{\tau}$ (see Figure 11a). Each vertical line connects the open and closed symbols for the same $Re_{\tau}$, indicating the difference between the corresponding $A_{1}$ and $C_{1}$ for each case within the groups. In Groups 1, 3, and 4, the estimated $C_{1}$ values differ substantially from their corresponding $A_{1}$ values. Groups 2 and 5 exhibit a spread of $A_{1}$ and $C_{1}$ values across different $Re_{\tau}$. Interestingly, while there is some variation in these values at lower $Re_{\tau}$, they appear to converge at higher $Re_{\tau}$, which seems to suggest that both $A_{1}$ and $C_{1}$ may depend on $Re_{\tau}$ and other atmospheric factors rather than being constants within the ASL. Figure 11b further explore the relation between $A_{1}$ and $C_{1}$, showing a weak-to-no correlation between the fitted $A_{1}$ and the estimated $C_{1}$ values among all the groups. Therefore, the data in this work do not support the $C_{1}=A_{1}$ relation proposed by AEM.