The atmospheric boundary layer (ABL) height plays a key role in many atmospheric processes as one of the dominant flow length scales. However, a systematic quantification of the ABL height over the entire range of scales (i.e., with periods ranging from one minute to one year) is still lacking in literature. In this work, the ABL height is quantified based on highresolution measurements collected by a scanning pulsed Doppler LiDAR during the recent American WAKE experimeNt (AWAKEN) campaign. The high availability of ABL height estimates ( $\u22482200$ collected over one year and each of them based on 10min averaged statistics) allows to robustly assess five different ABL height models, i.e., one for convective thermal conditions and four for stable conditions. Thermal condition is quantified by a stability parameter spanning three orders of magnitude and probed by nearground 3D sonic anemometry. The freeatmosphere stability, quantified by the Brunt–Väisälä frequency, is both calculated from simultaneous radiosonde measurements and obtained from the best fit of two of the chosen ABL height models. Good agreement is found between the data and three of the chosen models, quantified by mean absolute errors on the ABL height between 281 and 585 m. Furthermore, the seasonal variability of the convective ABL height model parameters (−15% to $+23%$ with respect to the year baseline) agrees with the variability of buoyancygenerated turbulence caused by the variation in solar radiation throughout the year.
I. INTRODUCTION
The atmospheric boundary layer (ABL), i.e., the lowest portion of the Earth's atmosphere affected by the presence of the ground via surface drag and energy exchange,^{1,2} is characterized by complex turbulence processes impacting several fields such as pollutant transport,^{3} numerical weather prediction,^{4} air quality,^{5} and wind^{6} and solar energy.^{7} A deep understanding of the physical processes within the ABL is rooted in reliable experimental measurements up to the ABL height [ $zi\u2248O(1000)$ m], which represents one of the fundamental scaling parameters governing the turbulence dynamics in the atmosphere.^{1,8,9} A key feature of this parameter is its large variability over a wide range of timescales, spanning from minuteslong^{1} to several months^{10} and years.^{11} To address this feature, accurate predictive models of the ABL height are available in the literature to capture the effect of either positive^{12,13} or negative^{14–18} vertical heat flux, as well as freeflow stratification in the troposphere,^{19,20} highlatitudes,^{21} and complex terrain effects. A caveat of these studies, however, is the limited amount of highquality experimental data (on the order of few tens of data points along the vertical direction,^{22} typically collected by 3D sonic anemometry) against which the ABL height models are assessed, mainly due to the relatively short duration of the experimental campaigns (e.g., a few months).^{23,24}
From the experimental standpoint, a suitable way to determine the ABL height involves diagnosing vertical turbulent statistics evolving along the vertical direction (z) and rapidly changing across the ABL height, and calibrating some empirical models dependent on z and using z_{i} as a free parameter.^{25} Several diagnostic methods are available in the literature according to the physical mechanisms driving the ABL evolution throughout the day. In particular, the daytime ABL height can be determined based on mean thermodynamic profiles,^{26} mean backscatter coefficient (which keeps track of the local aerosol concentration),^{27–31} backscattertoextinction coefficient ratio,^{32} backscatter coefficient variance,^{33} vertical velocity variance,^{34} and turbulent shear stress,^{35} i.e., statistical quantities reflecting a high turbulence intensity within the ABL and a subsequent sharp decrease moving into the free troposphere. By contrast, the nocturnal ABL, constituted by a thin shear layer (of the order of few hundreds of meters) closer to the ground and capped by a residual layer aloft, is more challenging to probe due to the low turbulence level and the presence of noncanonical ABL events. Thus, retrieval methods for the stable ABL height seek local maxima of the mean velocity^{36,37} as well as local minima of streamwise^{37} or vertical velocity^{38} fluctuations.
In this scenario, a breadth of experimental instruments are commonly used to determine the ABL height experimentally, such as radiosondes,^{39,40} ceilometers,^{31} and optical sensing techniques.^{10,26,27,40–42} Among the latter, the scanning pulsed Doppler light detection ranging (LiDAR) technology has gained importance in estimating the ABL height due to its extended measurement range (over 2000 m) together with its high spatiotemporal resolution ( $\u224820$ m, $\u22481$s, respectively) (see, e.g., Refs. 25, 30, and 43–46). Several diagnostic techniques have been developed to quantify the ABL height based on LiDARretrieved quantities, such as wavelet covariance transform,^{10,28–30} image processing,^{47} empirical models,^{27,48} and gradient^{44} and variance methods.^{33} However, the accuracy of these methods relies upon months to yearslong experimental campaigns to quantify the variability of z_{i} with adequate statistical convergence.
To address several research questions in the wind energy field, the multiinstitutional American WAKE experimeNt (AWAKEN) campaign was recently conducted in the Southern Great Plains (SGP) in Oklahoma, U.S., from 2022 to 2024.^{49,50} One of the main goals of this campaign is to study the mutual interaction between the ABL and five closely spaced wind farms under different micrometeorological conditions. This goal, in turn, requires a systematic and reliable quantification of the ABL height over the entire range of energycontaining turbulent scales. Thus, in this work, a novel dataset of ABL height measurements is presented based on data collected by one groundbased scanning LiDAR deployed during the AWAKEN campaign. In particular, highfrequency, qualitycontrolled LiDAR measurements of vertical velocity and signaltonoise ratio are used to estimate z_{i} under a wide range of thermal and shear conditions, the latter quantified by highfrequency singlepoint 3D sonic anemometer measurements close to the ground. Given the large number of available ABL height estimates, the present analysis aims to obtain a robust assessment of five ABL height predictive models available from the literature (both in stable and convective conditions) against the experimental ABL height measurements. Two of these models encompass the effect of the freeatmosphere stability, quantified by the Brunt–Väisälä frequency. The latter is then validated against the analogous estimates obtained from simultaneous and colocated radiosonde profiles of potential temperature. Finally, leveraging the yearlong span of the present measurements, the variability of the ABL height model parameters is linked to the presence of synoptic timescales driving the flow throughout the year.
This paper is organized as follows: in Sec. II, the experimental campaign, LiDAR scanning strategy, and quality control processes are detailed. In Sec. III, both the predictive and the diagnostic ABL height models used in this work are introduced. Results are presented in Sec. IV, and conclusions are discussed in Sec. V.
II. EXPERIMENTAL DATA COLLECTION
A. Experimental site, instrumentation, and scanning strategy
The data presented in this work are collected at the AWAKEN C1a site located midway between the King Plains and Armadillo Flats wind farms (cf. with Fig. 1), mainly featuring northsouth wind conditions.^{11,50–52} Among the instruments deployed at site C1a by different research teams, this study focuses on the data collected by two sensing instruments, namely, one Streamline XR scanning pulsed Doppler LiDAR (manufactured by Halo Photonics) and one surface flux station, both owned by The University of Texas at Dallas and deployed from October 4, 2022, to October 14, 2023. The surface flux station is equipped with two CSAT33D sonic anemometers (manufactured by Campbell Scientific Inc.) mounted at 2 m height along the prevailing wind direction (north–south) probing the threedimensional wind vector and temperature within a measurement volume of the order of 0.1 m and sampling rate of 20 Hz, thus fully resolving the shear and buoyancy turbulent motions close to the ground.
In order to probe vertical distributions of velocity variance ( $w\u2032w\u2032\xaf$) and timeaveraged signaltonoise ratio (SNR), which represents two quantities typically utilized to estimate the ABL height,^{32,53} the Streamline XR LiDAR performed 10minlong verticalstaring measurements with 48 m range gate and 2 Hz sampling rate. The measurement range spans from 120 m above the ground (end of the LiDAR's blind zone) up to 5976 m through 123 nonoverlapping gates. In total, 6734 vertical stare measurements (each of them with a 10min measurement period) were collected. For further details about the LiDAR's working principle, as well as the internal laser architecture, the reader is referred to the study of Wu et al.^{54}
As detailed in Sec. III, one of the fundamental scaling parameters governing the ABL evolution is the mean potential temperature gradient in the free atmosphere,^{2} which is quantified by the Brunt–Väisälä frequency. The latter will be estimated from the ABL height models calibrated onto the LiDAR data (see Sec. IV B). To validate the bestfitted values of the Brunt–Väisälä frequency, the radiosonde data collected by the DOE Atmospheric Radiation Measurement (ARM) at the SGP central facility (site C1, nearly 30 km north of site C1a, reported in Fig. 1 with white square symbol)^{55} are included in the present analysis. Each radiosonde returns instantaneous profiles of dry air temperature and pressure sampled with variable vertical resolution (0.1–40 m), thereby allowing for the calculation of the instantaneous profile of potential temperature from the ground proximity ( $\u2248300$ m) to well above the ABL height ( $\u22655000$ m). Five different hours across the day are considered in this study when radiosonde data are available: 02:00 UTC (21:00 local time), 05:00 UTC (00:00 local time), 11:00 UTC (06:00 local time), 17:00 UTC (12:00 local time), and 23:00 UTC (18:00 local time). Each launch is repeated across different days within the period October 4, 2022–October 14, 2023, in order to provide statistical significance to the retrieved quantities (50, 171, 188, 160, and 208 launches are done at each mentioned hour, respectively). For each considered time, the instantaneous vertical profiles of potential temperature are ensemble averaged over a common vertical coordinate spanning from 300 to 3000 m with vertical resolution of 50 m, and eventually utilized to calculate the vertical profiles of the Brunt–Väisälä frequency [see Eq. (5)].
B. Quality control of the sonic data
To obtain reliable estimates of these two statistical quantities, the following requirements are enforced:

The statistical uncertainty of $u\tau $ ( $\Delta u\tau $) must be lower than 20% of the expected value [calculated via Eq. (1)].

Krishnamurthy et al.^{11} have shown that hourly averaged values of friction velocity collected at the ARM SGP C1 site are statistically bounded within $0.25$ and $0.45\u2009 m s\u22121$ throughout the day. Thus, only friction velocity values smaller than $1 m s\u22121$ are considered in the rest of this study.

The statistical uncertainty of L ( $\Delta L$) must be lower than 50% of the expected value [calculated via Eq. (2)].
The present quality control of 3D sonic anemometer data is assessed against the Monin–Obukhov similarity theory; furthermore, details are reported in Appendix B. It is noteworthy that the statistical uncertainties of both $u\tau $ and L are calculated via bootstrap algorithm over 100 subset periods containing 500 points each^{57} using Eqs. (1) and (2), respectively. As expected, the enforcement of these requirements reduces the data availability. The relative percentage of data reduction introduced by each criterion is reported in Fig. 2(a) with black bars, while the cumulative percentage of data rejection is reported in Fig. 2(b). Here, it is observed that the largest source of data rejection (34.2%) is given by the unavailability of the sonic anemometer data, while the aforementioned criteria cause relatively small percentage of data reduction (9.9% and 13%, respectively). Eventually, of the 6734 LiDAR and sonic datasets, 2892 (42.9%) are left for further analysis after the quality control of the sonic anemometer data.
C. Quality control of LiDAR data
Focusing on the Streamline XR datasets, only LiDAR data taken simultaneously with the qualitycontrolled sonic anemometer data are considered for further analysis. Among them, those featuring cloudy conditions are removed. This requirement is checked based on the variance of the attenuated backscatter ( $\beta att$) at each height [ $\sigma \beta (z)$]. In particular, datasets showing $\sigma \beta \u226510\u22128m\u22122 s rad\u22122$ at any measured height are considered to be affected by the presence of clouds are discarded from further analysis. The attenuated backscatter is automatically calculated by the Streamline XR LiDAR model as the product of a calibration curve [ $C\beta (z)$] and the rangecorrected SNR ( $z2SNR$):^{58} $\beta att(z,t)=C\beta (z)z2SNR(z,t)$. Of the 2892 selected LiDAR datasets, 2274 (representing the 33.8% of the 6734 datasets initially available) are retained after the sonic and LiDAR quality control process (cf. with Fig. 2), which are subsequently utilized to estimate the ABL height.
For all the cloudfree LiDAR datasets, the dynamic filtering algorithm developed by Beck and Kühn^{59} is implemented to remove instantaneous outliers present in the vertical velocity time series, which could bias the quantification of this statistical moment. Specifically, for each time record at each height, instantaneous samples featuring either vertical velocity fluctuations outside of the interval ±3 m/s or SNR outside of $[\u221220,\u20090]$ dB are marked as outliers and removed. As shown in Sec. III B, the velocity variance profiles are utilized to estimate the ABL height only during nighttime hours, i.e., when the thermal conditions are mostly stable. Thus, although the choice of limiting the vertical velocity fluctuations to ±3 m/s might seem too restrictive for daytime convective conditions, this threshold does not influence the quantification of ABL height during nighttime, when the vertical velocity fluctuations are typically much less than ±3 m/s. Subsequent to the choice of thresholds, for each instantaneous realization, a normalized bivariate probability density function of w and SNR is calculated considering all the samples within 150 s. All the normalized $(w,SNR)$ samples with occurrence less than 0.1% are considered outliers, and the respective instantaneous velocity record is removed from the time series. Notably, this algorithm is implemented only to remove the incorrect estimates of vertical velocity, whereas no action is done on the SNR.
III. ATMOSPHERIC BOUNDARY LAYER HEIGHT QUANTIFICATION
A. Predictive models for the atmospheric boundary layer height
As mentioned in the Introduction, several models have been developed to predict the ABL height involving both nearground shear and thermal conditions and the free troposphere stability. In this subsection, five ABL height models available in the literature are introduced (one related to convective thermal conditions and four related to stable conditions), which will be subsequently validated against the ABL height estimates (Sec. IV B). The advantages and limitations of each of these models are discussed, as is the physical meaning of the free parameters utilized by the models.
1. The Rossby and Montgomery model
2. The Zilitinkevich and Mironov model
3. The Steeneveld model
The constant k is empirically quantified by Steeneveld et al.,^{18} plotting the lefthand side of Eq. (7) vs the argument in parenthesis for different intervals of N/f. For the present study, it is assumed that N/f in stable conditions is of the same order of magnitude as in convective conditions ( $N/f\u2248100$ at the SGP sites^{62,63}). The constant k is iteratively chosen based on this assumption; eventually, a value of $k=4\xd710\u22123$ is set in the remainder of this work.
It is noteworthy that all the bestfitted values mentioned in this work have been obtained in Matlab through the nonlinear least squares optimization method implementing a TrustRegion algorithm.^{64} The only constraint imposed on the optimization problem is that all the ABL model parameters must be greater than 0. The minimum of the cost function is then sought iteratively through the minimization of the secondorder local approximation of the cost function, where the Jacobian and Hessian matrices are evaluated trough secondorder central finite difference schemes. The iterative search is stopped when one of the following criteria is verified: a maximum of 400 iterations is reached; the cost function changes by a value less than $10\u22126$; the norm of the incremental step is less than $10\u22126$. Notably, for all the tested cases, the iterative algorithm stopped because the cost function increment was lower than $10\u22126$.
4. The Joffre model
5. The Kitaigorodskii model
B. Algorithms to estimate the atmospheric boundary layer height based on LiDAR data
In the context of the daytime convective ABL, the vertically staring Streamline XR Doppler LiDAR is utilized during the AWAKEN campaign to probe vertical turbulent motion. The timevsheight quantification of the attenuated backscatter ( $\beta att$), assumed as an indirect estimate of the local aerosol concentration,^{66} is provided by the LiDAR as well. The latter is a compelling quantity when it comes to estimating the convective ABL height, since the daytime turbulent mixing induces an almost uniform aerosol concentration throughout the mixed layer, whereas the entrainment zone aloft is characterized by a sharp decrease, eventually leading to the low aerosol concentration in the free atmosphere. Thus, the height associated with a rapid decrease in aerosol concentration (or, equivalently, $\beta att$) represents a physically sound estimate of the ABL height in the presence of a wellestablished mixed layer.^{32,67} It is noteworthy that, for the Streamline XR LiDAR used in the present work, the attenuated backscatter is calculated as^{58}^{,} $\beta att(z,t)=C\beta (z)z2SNR(z,t)$, where $C\beta (z)>0$ is a calibration curve provided by the manufacturer and accounting for the laser focus function. Since the latter may drift from the calibration conditions in time, a recalibration of $C\beta (z)$ against simultaneous and independent measurements of $\beta att$ is required. Lacking this possibility, for the present experiment, we assume the rangecorrected timeaveraged SNR ( $z2SNR$) as an indicator of the local aerosol properties (as previously done, e.g., by Kong and Yi,^{32} Menut et al.,^{33} and Toledo et al.^{68}).
Several diagnostic methods have been defined to capture rapid changes in the aerosol properties, for instance, based on the vertical evolution of the mean attenuated backscatter^{27} as well as the maximum attenuated backscatter variance,^{33} particle depolarization, and LiDAR ratio.^{32} Regardless of the particular quantity used, one method largely applied to individual sharp transitions of the aerosol properties is the wavelet covariance transform^{28,30,42,67,69} (WCT), which is used in this work to estimate the ABL height from LiDAR data during daytime. The details of this method are provided in Appendix A.
The WCT method, however, has two main limitations. First, it may lead to incorrect ABL estimates in the presence of multilayered aerosol structures. Second, a mixed layer must be well established in the ABL to maintain high and uniform aerosol concentrations. This is not the case during nighttime and early morning conditions, when the ABL is characterized by a thin layer of sheargenerated turbulence close to the ground and a residual layer aloft. To address the first point, a careful inspection of all the datasets is conducted based on the visualization of the timevsheight plots of vertical velocity and rangecorrected SNR; examples of these plots are reported in Figs. 3(a), 3(b), 3(d), 3(e), 3(g), and 3(h) for qualitycontrolled datasets. In regard to the second point, the WCT method is used only with datasets collected after 16:00 UTC time (corresponding to 11:00 local time), where a convective mixing is assumed to be present; this choice of threshold will be further discussed in Sec. IV A.
When the dataset does not show a sufficiently high positive buoyancy production (typically before 16:00 UTC), the vertical distribution of rangecorrected SNR is not a reliable source to estimate the ABL height. However, shearinduced turbulence is still expected to be generated at the ground and rapidly decrease moving upward. In the absence of noncanonical stable ABL events (such as gravity waves, bores, and shear instabilities), literature results have identified the ABL height in several ways, such as the height at which the streamwise velocity variance reaches a minimum^{37} as well as the height at which the variance of vertical velocity is less than a certain percentage of its value at the ground.^{22} In this work, we assume the local minimum of vertical velocity variance (after the quality control procedure described in Sec. II C) as the best estimator of the ABL height in stable conditions. An example of this method is visualized in Fig. 3(f).
After estimating the ABL height for all the available datasets, a further quality check is implemented to avoid the presence of nonphysical estimates in the database. These outliers are associated with extreme over or underestimation of the ABL height; thus, the z_{i} estimates are grouped into nonoverlapping hourlong bins and, for each of them, datasets featuring ABL height less than the 5th percentile or greater than the 95th percentile are further investigated. In particular, the ABL height is manually overwritten based on the visual inspection of timevsheight plots of vertical velocity and rangecorrected SNR. An example of this procedure is reported in Figs. 3(g)–3(i), where two layers of high aerosol concentration (roughly located at 1400 and 2200 m, respectively) are evidenced in the rangecorrected SNR pattern [Fig. 3(h)]. The WCT initially quantifies the ABL height as 2426 m (dashed line); however, the visual inspection of the vertical velocity pattern [Fig. 3(g)] shows negligible turbulent fluctuations above 1357 m (continuous line). Thus, the latter is assumed as the new ABL height value for the considered dataset. In total, 155 ABL height values (6.8%) are manually checked.
IV. RESULTS
A. Assessment of the daily cycle of atmospheric boundary layer height
After estimating the ABL height through the abovementioned diagnostic methods, the calculated z_{i} values are used to evaluate the daily cycle of ABL height. The latter is reported in Fig. 4(a) overlapped to the daily cycle of vertical velocity variance. The analogy between the distributions of z_{i} and $w\u2032w\u2032\xaf$ confirms that the primary mechanism governing the ABL evolution at the tested site is the thermal stability. From Fig. 4(a), the vertical velocity variance reaches its maximum between 16:00 and midnight (UTC time) due to the onset of a positive buoyancy flux which, in turn, enhances the aerosol mixing within the ABL. As anticipated in Sec. III B, the WCT method can only be applied in the presence of significant turbulent mixing; thus, it is utilized only for datasets collected between 16:00 and midnight, whereas the minimum vertical velocity variance method is used for the remaining data. This time cutoff is highlighted in Fig. 4(a) with a vertical dashed line.
The ABL height can be further assessed comparing the current daily cycle against the analogous estimates obtained by Krishnamurthy et al.,^{11} who used a machine learning algorithm on a yearslong scanning LiDAR database collected at the ARM SGP C1 site ( $\u224830$ km north of the C1a site) via verticalstaring scans. The result, reported in Fig. 4(b), demonstrates strong quantitative agreement between the current daily cycle and the literature result, both in the hourly averaged ABL height values (black symbols) and in the standard deviation (vertical uncertainty interval). The temporal oscillation of the current daily cycle is deemed to be caused by the relatively short sampling period (10 min) of each vertical scan.
B. Stability dependence of the atmospheric boundary layer height
After comparing the current ABL height estimates against previous results, the dependence of z_{i} on the thermal conditions is investigated by plotting the stability parameter $zi/L$ against μ (Fig. 5) distinguishing between convective [Fig. 5(a)] and stable [Fig. 5(b)] conditions. To facilitate the visual comparison between the data and the models, as well as the departure of the data from the $1:1$ linear trend [dashed lines in Figs. 5(a) and 5(b)], binaveraged values of $zi/L$ (white filled squares) are overlapped on the data.
1. Convective conditions
Focusing on convective conditions first [Fig. 5(a)], strong agreement is found between the data and the calibrated model of Kitaigorodskii and Joffre^{12} [Eq. (10)], the latter fitted over 1078 data points. The relative error distribution between the binaveraged data values and the Kitaigorodskii and Joffre^{12} model prediction, reported in Fig. 5(c) with red symbols and continuous line, indicates percentage errors between −19% and 15%. However, the rootmeansquare error (RMSE) and mean absolute error (MAE) between the ABL height evaluated experimentally and predicted by the model are quite large (798 and 585 m, respectively, as reported in Table I); these large values are deemed to be caused both by the large variability of the convective ABL height values [visualized by the large scattering of data points in Fig. 5(a) and by the large uncertainty intervals in Fig. 4(b)] and by the lack of direct quantification of the Brunt–Väisälä frequency simultaneous to each ABL height estimate.
Reference .  Eq. # .  Fit parameters .  z_{i} RMSE (m) .  z_{i} MAE (m) . 

Ref. 12  10  $m\u2033=9.31\u2009(\xb11.42)$  798  585 
$N/f=90\u2009(\xb18)\u2009(L<0)$  
Ref. 18  7  $\alpha =3.0\u2009(\xb10.2)$  474  281 
$N/f=85\u2009(\xb113)\u2009(L>0)$  
Proposed model  9  $C=0.119\u2009(\xb10.017)$  484  288 
$\beta =2.7\xd710\u22123\u2009(\xb14.6\xd710\u22123)$  
$\gamma =1.22\u2009(\xb10.28)$ 
Reference .  Eq. # .  Fit parameters .  z_{i} RMSE (m) .  z_{i} MAE (m) . 

Ref. 12  10  $m\u2033=9.31\u2009(\xb11.42)$  798  585 
$N/f=90\u2009(\xb18)\u2009(L<0)$  
Ref. 18  7  $\alpha =3.0\u2009(\xb10.2)$  474  281 
$N/f=85\u2009(\xb113)\u2009(L>0)$  
Proposed model  9  $C=0.119\u2009(\xb10.017)$  484  288 
$\beta =2.7\xd710\u22123\u2009(\xb14.6\xd710\u22123)$  
$\gamma =1.22\u2009(\xb10.28)$ 
The calibrated model parameters ( $m\u2033,\u2009N/f$) are reported in Table I as well. The fitted value of $m\u2033\u2009(9.31)$ is of the same order of magnitude as the one proposed by Kitaigorodskii and Joffre^{12} ( $m\u2033=6$), although larger. A possible justification of this discrepancy may be found in the synoptic variability of the ABL height across different seasons,^{11} which will be investigated in Sec. IV C. Finally, the calibrated value of $N/f=90$ agrees with previous estimates of the nondimensional Brunt–Väisälä frequency ( $N/f\u2248102$) at SGP latitudes.^{62} Furthermore, the vertical profiles of the Brunt–Väisälä frequency obtained from the radiosonde data from site C1 are reported in Fig. 6(f) for each considered hour when launches are available. The vertical profiles of mean potential temperature utilized to calculate N/f are reported in Figs. 6(a)–6(e) for each considered launch time. The N/f interval evaluated from radiosondes during daytime convective conditions [launches at 17:00 and 23:00 UTC, reported with red and magenta lines in Fig. 6(f)] spans from 50 to 200 at different heights, thus in good agreement with the value retrieved from the Kitaigorodskii and Joffre^{12} model ( $N/f=90$). In conclusion, the current estimates of the convective ABL height from AWAKEN data lead to an overall good level of agreement with the theoretical prediction of Kitaigorodskii and Joffre,^{12} although this model is not able to completely capture the ABL height variability in convective conditions.
2. Stable conditions
For stable conditions [Fig. 5(b)], we first observe an increase in $zi/L$ with slope greater than 1 increasing the stability ( $\mu \u2a86102$), consistently with previous results.^{13,18} Thus, the model of Zilitinkevich and Mironov^{16} [reported in Fig. 5(b) with blue line], which predicts a lessthanlinear power law for very stable conditions, largely underestimates the ABL height. Notably, this outcome does not depend on the uncertainty on the Brunt–Väisälä frequency as evidenced from the blue shaded area in Fig. 5(b), which is obtained solving Eq. (6) letting $1\u2264N/f\u2264103$. Similarly, the parameterization of Joffre et al.^{13} [Eq. (8), reported in Fig. 5(b) with cyan line), although it models the monotonic increase in $zi/L$ with μ, underestimates the stabilitycorrected ABL height by 26% or more as evidenced by the percentage error distribution in Fig. 5(c) (cyan symbols).
Better agreement is found calibrating the model of Steeneveld et al.^{18} [Eq. (7), green line in Fig. 5(b)], which leads to large percentage errors (−69% to $+51%$) for neutral to moderately stable conditions ( $\mu \u226450$) followed by smaller values ( $<10%$) for the rest of the stability interval [cf. with Fig. 5(c)]. This feature indicates that the model of Steeneveld et al.^{18} is better suited for strongly stable conditions. The overall fitted value of α (3.0, cf. with Table I), however, is in excellent agreement with the one proposed by Steeneveld et al.^{18} (α = 3), and the obtained value of N/f (85 ± 13) is in good agreement with previous estimates in the SGP area during nighttime conditions (N = 0.01 Hz, corresponding to $N/f\u2248116$).^{63} Similarly, the N/f interval in stable conditions retrieved from the radiosonde data ( $100\u2264N/f\u2264300$ indicated by the cyan line in Fig. 6(f), related to radiosonde launches done at 02:00 UTC) agrees well with the value calibrated from the model of Steeneveld et al.^{18}
Finally, the currently proposed parameterization of the stable ABL height [Eq. (9)] is calibrated onto the experimental data and reported in Fig. 5(b) with the magenta line. It is noticed that, among the models introduced to parameterize the stable ABL height, Eq. (9) leads to the smallest percentage error distribution across μ, mainly due to the fact that it relies on the largest number of free parameters. The value of the constant C (0.119, cf. with Table I), quantifying the linear growth rate of the ABL for weak stability, is in good agreement with the order of magnitude of C in Eq. (3) for stable conditions ( $C\u22480.12$) found in the literature.^{70} This confirms the initial hypothesis that Eq. (9) may be seen as a generalization of the geostrophic drag law valid for the entire stability range. The fitted value of $\beta =(2.7\xb14.6)\xd710\u22123$ is uncertain due to the large scattering of ABL height values across the stability range. This is related to the difficulty of estimating the ABL height while lacking a welldefined boundary between turbulent and nonturbulent portions of the atmosphere for nighttime conditions. Finally, the value $\gamma =1.22$ captures the growth of $zi/L$ for strongly stable conditions. It is noteworthy that Eq. (9) can be approximated by $zi/L\u2248C\beta \mu \gamma +1$ as $\mu \u226b1$. Given the current value of γ, we obtain the following expression for strongly stable conditions: $zi/L\u221d\mu 2.22$. Joffre et al.^{13} estimated a growth rate proportional to $\mu 1.5$ for $10\u2264\mu \u2264100$. The faster growth of $zi/L$ currently found is thought to be a consequence of the larger extension of the current stability range ( $1\u2264\mu \u2264103$) over which the behavior of $zi/L$ is modeled.
For stable conditions, the RMSE (484 and 474 m) and MAE (288 and 281 m) values between the experimental ABL height values and the model predictions are large (see Table I), although lower than what was found for convective conditions. It is noticed that the models of Eqs. (7) and (9) lead to similar values of RMSE and MAE. This feature, together with the agreement between each model's constants and the literature values, leads to the conclusion that the discrepancy between the experimental values and the prediction is due to the uncertainties in the experimental quantification of the ABL height rather than a limitation of the models to capture the flow physics. Thus, the current stable ABL height dataset is in good agreement with the stability trend predicted by the models, even though the uncertainty associated with the experimental ABL height estimates leads to significant mean error values with respect to the models' predictions.
C. Seasonal variability of the parameters
The availability of ABL height estimates over one year is now leveraged to quantify the effect of seasonal timescales on the ABL height model parameters. To this aim, the data are binned into four subsets spanning three months each, namely, December–January–February (DJF), March–April–May (MAM), June–July–August (JJA), and September–October–November (SON). The daily cycles of z_{i} and vertical velocity variance are reported in Fig. 7 for each data bin. Focusing on daytime convective conditions (i.e., after 16:00 UTC, reported with vertical dashed line), smaller velocity variance is observed for the period DJF with respect to the baseline daily cycle [Fig. 7(a)] due to the shorter exposure of the ground to the solar irradiation; this results in less buoyancygenerated TKE and thus lower ABL height values [blue symbols in Fig. 8(a)] with respect to the baseline distribution during daytime conditions. By the same principle, higher values of $w\u2032w\u2032\xaf$ (associated with higher values of ABL height) are found during MAM and JJA daytime periods [Figs. 8b and 8(c)], whereas negligible difference is found during SON. By contrast, before sunrise (i.e., before 16:00 UTC) no qualitative difference is found in the daily averaged values of z_{i} within the seasons with respect to the baseline trend.
The impact of seasonal trends on the ABL height is further examined by calibrating the models of Eqs. (7) and (9) for stable conditions and the model of Eq. (10) for convective conditions over different seasonal subsets. The fitted parameters are reported in Table II grouped by model and in Fig. 8 by physical effect. In particular, Fig. 8(a) depicts all the parameters directly impacted by any vertical displacement of the μ $zi/L$ distribution (namely, $m\u2033,\u2009\alpha $, and C), which is deemed to be caused by the seasonal timescale. Figure 8(b) reports the averaged values of N/f across different seasons for both stable and convective conditions. Finally, Fig. 8(c) focuses on the parameters $\beta ,\u2009\gamma $ of Eq. (9).
.  Ref. 12 (Eq. 10) .  Ref. 18 (Eq. 7) .  Present model [Eq. (9)] .  

.  $m\u2033$ .  $Nf$ (L < 0) .  α .  $Nf$ (L > 0) .  C .  β ( $\xd710\u22123$) .  γ . 
DJF  8.3 (–11%)  88 (–2%)  3.49 ( $+16%$)  64 (–24%)  0.13 ( $+9%$)  2.15 (–22%)  1.17 (–4%) 
MAM  11.5 ( $+23%$)  92 ( $+2%$)  2.32 (–23%)  80 (–6%)  0.14 ( $+22%$)  1.02 (–63%)  1.37 ( $+12%$) 
JJA  11.4 ( $+13%$)  100 ( $+11%$)  2.91 (–3%)  97 ( $+15%$)  0.11 (–6%)  1.34 (–51%)  1.45 ( $+19%$) 
SON  7.9 (–15%)  85 (–6%)  2.93 (–3%)  98 ( $+16%$)  0.12 (–1%)  1.31 (–52%)  1.41 ( $+15%$) 
.  Ref. 12 (Eq. 10) .  Ref. 18 (Eq. 7) .  Present model [Eq. (9)] .  

.  $m\u2033$ .  $Nf$ (L < 0) .  α .  $Nf$ (L > 0) .  C .  β ( $\xd710\u22123$) .  γ . 
DJF  8.3 (–11%)  88 (–2%)  3.49 ( $+16%$)  64 (–24%)  0.13 ( $+9%$)  2.15 (–22%)  1.17 (–4%) 
MAM  11.5 ( $+23%$)  92 ( $+2%$)  2.32 (–23%)  80 (–6%)  0.14 ( $+22%$)  1.02 (–63%)  1.37 ( $+12%$) 
JJA  11.4 ( $+13%$)  100 ( $+11%$)  2.91 (–3%)  97 ( $+15%$)  0.11 (–6%)  1.34 (–51%)  1.45 ( $+19%$) 
SON  7.9 (–15%)  85 (–6%)  2.93 (–3%)  98 ( $+16%$)  0.12 (–1%)  1.31 (–52%)  1.41 ( $+15%$) 
Focusing on Fig. 8(a), the parameter α of Eq. (7) (reported with squared symbols) assumes higher value on DJF with respect to the reference ( $+16%$, cf. with Table II), indicating a lower ABL height under stable conditions as compared to the averaged value across the year. By contrast, a sensibly lower α value (−23% from baseline) is found during MAM, followed by slightly lower values during JJA ((−3%) and SON ((−3%), which implies higher stable ABL height during these months. However, it should be noticed that the best fit for the MAM bin relies on fewer data points as compared to the other bins, as reported by Fig. 8(d). This implies larger fitting uncertainty and thus less reliable values for the MAM period. These considerations are less evident when observing C, for which the uncertainty intervals associated with the fitting do not allow for firm conclusions of the seasonal trend [cf. with Fig. 8(a), triangular symbols]. Finally, the seasonal distribution of $m\u2033$ (thus focused on convective conditions) confirms the previously mentioned qualitative result, i.e., lower convective ABL height during DJF with respect to the year average modeled by a lower value of $m\u2033$ (–11%), and higher ABL heights (thus higher $m\u2033$) during MAM ( $+23%$ increase in $m\u2033$) and JJA ( $+13%$ increase) and, finally, lower ABL heights during SON (–15%). In conclusion, it is inferred that the effect of seasonal timescales on the ABL height is embedded into the model parameters α (for stable conditions) and $m\u2033$ (for convective conditions), whereas the parameter C is less sensitive to the seasonal variability.
The seasonal trends of nondimensional Brunt–Väisälä frequency obtained by calibrating Eq. (7) in stable conditions and Eq. (10) in convective conditions exhibit overlapping intervals ( $64\u2264N/f\u226497$ in stable conditions and $85\u2264N/f\u2264100$ in convective conditions). However, while the Brunt–Väisälä frequency is nearly constant across different bins for convective conditions, the analogous frequency in stable conditions shows a peak during JJA and SON [as reported in Fig. 8(b)]. The uncertainty intervals obtained for the MAM data [green symbols in Fig. 8(b)] are sensibly larger than the other estimates, again due to the lower number of data points available.
Finally, the seasonal trend of γ [reported in Fig. 8(c) with diamond symbols] reveals a slight increase from DJF (–4% with respect to the baseline) to SON ( $+15%$). However, the large uncertainty intervals do not allow for further considerations about the growth rate of z_{i} during stable conditions across different seasons. Similarly, the uncertainty intervals around β do not allow for definitive conclusions on the effect of seasonal scales.
V. CONCLUSION
In this work, the ABL height (z_{i}) is quantified during the yearlong AWAKEN campaign. Based on time series of vertical velocity and rangecorrected signaltonoise ratio probed by a vertically pointing scanning pulsed Doppler LiDAR, the ABL height is systematically evaluated for over 2200 qualitycontrolled 10min averaged measurements. Two different diagnostic models are used to quantify the ABL height under convective and stable conditions. One model performs best under daytime conditions (i.e., when a wellmixed layer is ingrained in the ABL), and the other model performs best during nighttime lowturbulent conditions. This approach provides the complete coverage of the ABL height evolution throughout the diurnal cycle. The average daily variation agrees well with previous literature results and with the daily cycle of vertical velocity variance.
The simultaneous availability of two colocated 3D sonic anemometers allows us to provide quantitative information on the shear and buoyancygenerated turbulence near the ground, thereby relating the ABL height variability to the nearground thermal stability and shear turbulence. Similarly, the freeatmosphere stability (expressed by the Brunt–Väisälä frequency, N) is quantified based on the radiosonde potential temperature profiles available from the neighboring ARM SGP C1 site. Five stabilitydependent ABL height models are introduced to capture the ABL height variability over a wide range of stability classes, the latter ranging within three orders of magnitude both for stable and convective conditions. One model (first proposed by Kitaigorodskii and Joffre^{12}) is used to assess the ABL height in convective conditions, while four models are introduced for stable conditions. Among the latter, three models are taken from the literature, namely, from the Zilitinkevich and Mironov,^{16} Joffre et al.,^{13} and Steeneveld et al.^{18} studies, while the last model is introduced here to generalize the geostrophic drag law originally proposed by Rossby and Montgomery.^{60} It is noteworthy that the Brunt–Väisälä frequency is fitted on the data as a free parameter for the models of Kitaigorodskii and Joffre^{12} and Steeneveld et al.,^{18} and it is subsequently validated against radiosonde measurements collected at the ARM SGP C1 site (located nearly 30 km north of the tested site). Good agreement is inferred for the Kitaigorodskii and Joffre,^{12} Steeneveld et al.,^{18} and proposed model with respect to the experimental ABL height estimates. By contrast, low agreement is obtained between the retrieved ABL heights against the Zilitinkevich and Mironov^{16} and Joffre et al.^{13} predictions. From a quantitative standpoint, the models' parameters match well with the analogous values found in the literature. Similarly, the fitted values of the Brunt–Väisälä frequency in convective ( $N/f=90$, where f is the Coriolis parameter) and stable ( $N/f=85$) conditions agree well with the analogous intervals obtained from radiosonde measurements ( $50\u2264N/f\u2264200$ and $100\u2264N/f\u2264350$ for convective and stable stratification, respectively).
Finally, the variability of each model's fitting parameters is physically linked to intraannual timescales by fitting the models over nonoverlapping seasonal subsets. In particular, the lower TKE buoyancy production (due to shorter periods of solar irradiance on the ground) during the winter season induces lower ABL height, as opposed to spring and summer seasons. This effect is carried on by each model's multiplicative constant, for which the intraannual distributions are found to agree with the seasonal cycle of TKE buoyancy production.
ACKNOWLEDGMENTS
We would like to acknowledge the essential contribution of the AWAKEN science team led by the National Renewable Energy Laboratory and its institutional and industrial partners. We also acknowledge the support of the technical staff for installing and maintaining the instruments and the landowners from Noble and Garfield counties in Oklahoma for granting access to their land. Data were obtained from, and this research was supported by the Atmospheric Radiation Measurement (ARM) user facility, a U.S. Department of Energy (DOE) Office of Science user facility managed by the Biological and Environmental Research program. This work was authored in part by the National Renewable Energy Laboratory, operated by Alliance for Sustainable Energy, LLC, for the U.S. Department of Energy (DOE) under Contract No. DEAC3608GO28308. Funding provided by the U.S. Department of Energy Office of Energy Efficiency and Renewable Energy Wind Energy Technologies Office. The views expressed in the article do not necessarily represent the views of the DOE or the U.S. Government. The U.S. Government retains and the publisher, by accepting the article for publication, acknowledges that the U.S. Government retains a nonexclusive, paidup, irrevocable, worldwide license to publish or reproduce the published form of this work, or allow others to do so, for U.S. Government purposes. LLNL is operated by Lawrence Livermore National Security, LLC, for the DOE, National Nuclear Security Administration under Contract DEAC5207NA27344. NREL is a national laboratory of the U.S. Department of Energy, Office of Energy Efficiency and Renewable Energy, operated by the Alliance for Sustainable Energy, LLC.
The authors are grateful to Raghavendra Krishnamurthy (PNNL) for providing the statistics of ABL height from the ARM database [Fig. 4(b)].
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
Author Contributions
M. Puccioni: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Writing – original draft (equal); Writing – review & editing (equal). C. F. Moss: Investigation (equal); Methodology (equal). M. S. Solari: Investigation (equal); Methodology (equal). S. Roy: Investigation (equal); Methodology (equal). G. V. Iungo: Funding acquisition (equal); Investigation (equal); Methodology (equal); Supervision (equal). S. Wharton: Funding acquisition (equal); Project administration (equal); Supervision (equal). P. Moriarty: Funding acquisition (equal); Project administration (equal); Supervision (equal).
DATA AVAILABILITY
The lidar data collected during the precampaign test are available at www.a2e.energy.gov.
APPENDIX A: THE WAVELET COVARIANCE TRANSFORM METHOD
In this work, b is assumed equal to the LiDAR range gate values within $zmin=500$ and $zmax=3500$ m, thus with increase equal to the range gate (48 m).
As evidenced by Eq. (A3), the correct quantification of z_{i} relies on a correct choice of dilation parameter, a, which is physically associated with the depth of the entrainment zone. In particular, a relatively low value of a ( $\u2248100$ m) is sensitive to the noise fluctuations of $z2SNR$, whereas a large value of a ( $\u2248103$ m) may exceed the entrainment zone depth, and thus, the local maximum of the WCT would be smoothed out by the convolution integral. In this work, the correct value of a is found by the iterative procedure introduced by Brooks.^{30} A large dilation parameter ( $a0=2000$ m) is initially guessed to obtain a first estimate of the WCT (said $WCT0$), the latter attaining a peak value ( $WCTmax,0)$ for a certain translation parameter, b_{0}, and reducing to a fraction of the maximum value (here assumed $0.3WCTmax,0$) within a certain interval $\Delta b0$ centered around b_{0} [cf. with Fig. 8(b)]. Since the WCT is nearly zero outside of the entrainment zone, $\Delta b0$ is assumed as a proxy of the entrainment zone depth^{30} and, thus, can be used as a better estimate of the dilation parameter for the next iteration ( $a1=\Delta b0$). The iterative procedure ends when the percentage difference between a_{i} and $\Delta bi\u22121$ (where i is the iteration number) is less than 10% or when the number of iterations reaches 100.
APPENDIX B: ASSESSMENT OF THE QUALITY CONTROL OF 3D SONIC DATA
The best fit of Eq. (B1) (whose parameters are reported in Table III for each velocity component) on the present qualitycontrolled data captures the increase in shearcorrected standard deviation with the stability parameters, as well as the plateau as z/L goes to 0. Comparing the bestfitted quantities with the reference values, the largest source of discrepancy is found in the values of b (Table III). This could be explained considering the lack of data beyond $z/L=2$, as opposed to the dataset shown by Pahlow et al. reaching $z/L=10$; this leads to a less precise calibration of the second term on the righthand side of Eq. (B1), which models the increase in shearcorrected standard deviation with the stability. However, the 3D sonic velocity statistics remaining from the quality control procedure agree with the MOST prediction, whereas the rejected portion of the data does not. Thus, it is concluded that the adopted quality control process leads to a reliable quantification of the turbulence state close to the ground.
Equation (B1) parameters .  a .  B .  c .  

Velocity component .  Ref. 71 .  Present .  Ref. 71 .  Present .  Ref. 71 .  Present . 
Streamwise (u)  2.3  2.1  4.3  14.4  0.5  0.6 
Spanwise (v)  2.0  1.0  4.0  13.6  0.6  0.5 
Vertical (w)  1.1  0.7  0.9  7.3  0.6  0.7 