Refining Planetary Boundary Layer Height Retrievals From Micropulse‐Lidar at Multiple ARM Sites Around the World

Knowledge of the planetary boundary layer height (PBLH) is crucial for various applications in atmospheric and environmental sciences. Lidar measurements are frequently used to monitor the evolution of the PBLH, providing more frequent observations than traditional radiosonde‐based methods. However, lidar‐ derived PBLH estimates have substantial uncertainties, contingent upon the retrieval algorithm used. In addressing this, we applied the Different Thermo‐Dynamic Stabilities (DTDS) algorithm to establish a PBLH data set at five separate Department of Energy's Atmospheric Radiation Measurement sites across the globe. Both the PBLH methodology and the products are subject to rigorous assessments in terms of their uncertainties and constraints, juxtaposing them with other products. The DTDS‐derived product consistently aligns with radiosonde PBLH estimates, with correlation coefficients exceeding 0.77 across all sites. This study delves into a detailed examination of the strengths and limitations of PBLH data sets with respect to both radiosonde‐ derived and other lidar‐based estimates of the PBLH by exploring their respective errors and uncertainties. It is found that varying techniques and definitions can lead to diverse PBLH retrievals due to the inherent intricacy and variability of the boundary layer. Our DTDS‐derived PBLH data set outperforms existing products derived from ceilometer data, offering a more precise representation of the PBLH. This extensive data set paves the way for advanced studies and an improved understanding of boundary‐layer dynamics, with valuable applications in weather forecasting, climate modeling, and environmental studies.


Introduction
The planetary boundary layer (PBL) is the lowest section of the troposphere, significantly influenced by terrestrial surfaces.This layer reacts to surface forcings, such as frictional drag and heat transfer, and is instrumental in regulating pollutant emission and dissipation (Caughey, 1984;Holtslag & Nieuwstadt, 1986;Z. Li et al., 2017;Mahrt, 1999;Monks et al., 2009;Stull, 1988).Interactions between the atmosphere and the surface culminate in substantial variations in the PBL Height (PBLH) across both spatial and temporal dimensions (Cimini et al., 2020;Guo et al., 2016;Illingworth et al., 2019;Kalmus et al., 2022;Seidel et al., 2010).Temporal fluctuations can emerge within an hour, while spatial changes may span several hundred meters to a few kilometers.The PBL exerts a pivotal influence on numerous atmospheric processes and applications, ranging from the dispersion of air pollutants to weather forecasting and climate modeling (Chu et al., 2019;Guo et al., 2021;Holtslag & Boville, 1993;Hu et al., 2010;Knote et al., 2015;Z. Li et al., 2017;Seidel et al., 2010;Su et al., 2024).Hence, the precise estimation of PBLH is vital to advancing our understanding of surface-atmospheric processes and refining climate and weather simulations.Yet, the measurement and modeling of PBLH pose significant challenges due to its high variability and limitations inherent in observational technology (Dang et al., 2019;Kotthaus et al., 2023; Traditionally, the estimation of PBLH has relied on in-situ radiosonde (RS) data, capturing thermodynamic profiles of temperature, humidity, and wind across various pressure levels (H.Li et al., 2021;Seidel et al., 2010;W. Zhang et al., 2018).Methods used in conjunction with RS data include surface-based inversions, relative humidity, and potential temperature gradient methods, Richardson number methods, and the parcel method (Bradley et al., 1993;S. Liu & Liang, 2010;Seidel et al., 2010;Vogelezang & Holtslag, 1996), among others.
While RSs have been employed as a standard in estimating PBLH, their precision varies notably depending on atmospheric thermodynamic conditions, with no single method proving universally effective across all conditions (H.Li et al., 2021).Furthermore, RS retrievals face restrictions in terms of spatial coverage and observational frequency, typically occurring only twice daily at most stations.In recent decades, lidar techniques have emerged as a potent alternative for PBLH estimation, offering a very high frequency, albeit with a lesser degree of accuracy compared to RS techniques (Eresmaa et al., 2006;B. Liu et al., 2018;Sawyer & Li, 2013;W. Zhang et al., 2016;D. Zhang et al., 2022).Lidar beams trace the diurnal evolution of the PBLH, measuring the backscattered radiation from aerosols confined within the PBL (Cohn & Angevine, 2000;Eresmaa et al., 2006;Hageli et al., 2000).Ideally, aerosols within the PBL create a pronounced contrast with the free atmosphere, facilitating the estimation of PBLH (Brooks, 2003;Campbell et al., 2002;Davis et al., 2000;Flamant et al., 1997).However, lidars may struggle to estimate PBLH in stable conditions or when there are multiple aerosol layers (Kotthaus et al., 2023;D. Zhang et al., 2022).
While the Different Thermo-Dynamic Stabilities (DTDS) method (Su et al., 2020) has showcased reasonably good accuracy in PBLH estimation at the Southern Great Plains (SGP) site, it is crucial to authenticate its performance in other regions.Moreover, the generation of a more precise PBLH data set for these areas can deepen our understanding of the PBL's behavior across diverse geographical regions with distinct climate regimes.As a result, this study aims to examine the DTDS algorithm's efficacy in estimating PBLH at five distinct Department of Energy (DOE) Atmospheric Radiation Measurement (ARM) sites.These comprise the SGP site and four ARM Mobile Facilities, namely the Green Ocean Amazon-GoAmazon (MAO); the Cloud, Aerosol, and Complex Terrain Interactions-CACTI (COR); the Convective and Orographically-Induced Precipitation Study (FKB), and the Biogenic Aerosols-Effects on Clouds and Climate-BAECC (TMP).This paper is organized as follows.Section 2 briefly elucidates the ground-based measurements at the ARM sites and the methodologies employed for calculating the PBLH.Section 3 presents the results, further subdivided into an overview of the PBLH product retrieved using the DTDS algorithm (hereinafter referred to as DTDS-PBLH), an assessment of DTDS-PBLH's performance against RS PBLH (RS-PBLH), comparison with another lidarbased PBL product using ceilometer data and investigate their errors and limitations.The final section provides a concise discussion and concludes the paper.

Data and Methodology
The PBLH varies significantly from site to site according to multiple climatic and meteorological conditions, both of which are highly determined by the geographical location of each site.To comprehensively evaluate the performance of the DTDS algorithm, the present study estimates the PBLH over five different DOE ARM sites, including one permanent location and four mobile facilities.ARM ground observation facilities have provided continuous and comprehensive measurements of all meteorological variables in the atmosphere and at the surface since the early 1990s.Specifically, this study draws information from the SGP site in the United States, the MAO site of the GoAmazon campaign in Brazil, the COR site of the CACTI campaign in Argentina, the FKB site of the COPS campaign in Germany, and the TMP site of the BAECC campaign in Finland.These sites were selected to address different scientific objectives, but their observational configurations are more or less standard, including the characterization of boundary-layer meteorology.Regarding surface characteristics, we only consider continental sites not influenced by a marine environment because the original DTDS algorithm was developed using data collected at the SGP site.Also, only those sites with more than 6 months of data are considered.Figure 1 shows the geographical distribution of the five sites selected.Table 1 summarizes each site's specifications, including their latitude, longitude, altitude, period of observations, and the cloud product used in this study.

Estimation of the PBLH Using RS Measurements
The ARM SONDE product includes atmospheric thermodynamic profiles obtained from its balloon-borne sounding system.These profiles contain height/pressure, temperature, humidity, wind speed, and direction data, captured at least four times daily (at 05:30,11:30,17:30,and 23:30 UTC) at each of the five ARM sites selected.Holdridge (2020) provide details about the SONDE product.This study used only RS data collected between 06:30 and 19:00 Local Time (LT).Because the vertical resolution of RS measurements fluctuates in accordance with the balloon's ascent rate, we resampled RS data through linear interpolation to attain a vertical resolution of five hPa.
We calculated the PBLH from radiosonde observations (i.e., RS-PBLH) using the well-established S. Liu and Liang (2010) method.The method starts by distinguishing the thermodynamics of the PBL by comparing the nearsurface potential temperature difference with a stability threshold δ s defined as 1 k for continental regions.Under convective conditions, this method identifies the PBLH as the height at which an air parcel rising adiabatically becomes neutrally buoyant.The method first scans upward a level k, where the temperature difference θ k θ 1 exceeds a threshold value δ u , which equals 0.5 k.This level is corrected by a second scan searching for the overshooting level, which defines the entrainment zone.The Liu and Liang method defines this overshooting level as the height where ∂θ k /∂z ≥ θr , with θr equal to 4 k/km.Under stable conditions, the PBLH is harder to determine because turbulence arises from buoyancy or wind shear.Liu and Liang's method defines a stable PBLH as either the top of the stable layer above the ground or at the low-level-jet nose height, whichever is lower.

Estimation of the PBLH Using Micropulse Lidar Measurements and the DTDS Algorithm
Operating at 532 nm in both parallel and perpendicular polarizations, the autonomous, ground-based micropulse lidar (MPL) provides backscattered radiation profiles at a temporal resolution of 10-30 s and a vertical resolution of 30 m. Raw data are corrected for background subtraction, signal saturation, overlap, after-pulse, and range per standard lidar-data processes (Campbell et al., 2002(Campbell et al., , 2003)).Quality-control flags help to omit erroneous data.
In principle, an inversion layer caps the PBLH, trapping moisture and aerosols within the boundary layer, thus creating a sharp contrast with the free atmosphere (Seibert, 2000).This contrast, visible in the backscattered radiation observed by lidar instruments, makes them ideal for estimating the PBLH.However, estimating the PBLH from lidar data is challenging when there is strong atmospheric stratification or when the inversion layer resides beneath the lidar blind zone.These situations occur more frequently during stable boundary-layer conditions (Z.Li et al., 2017).Compared to the previous approaches.DTDS methods have made several refinements.First, DTDS incorporates both the wavelength and gradient methods.In addition, DTDS considers the diurnal variations of the boundary layer, allowing for a more comprehensive estimation of the PBLH with enhanced vertical consistency and temporal continuity, distinguishing it from the traditional method that relies solely on a single backscatter profile and does not consider the diurnal cycle.
Traditional methods face challenges when estimating the PBLH under cloudy conditions, leading to potential cloud contamination.To address this issue, the DTDS approach effectively handles cloudy conditions by assessing the cloud-surface coupling.DTDS uses Lidar-derived PBLH, cloud position, and the lifted condensation level (LCL) to classify clouds as coupled or decoupled to the surface.This differentiation of clouds aids in  accurately diagnosing the PBLH (Su et al., 2022).All these refinements in DTDS offer more reliable and robust PBLH estimations, even in complex atmospheric conditions.
In this study, minor modifications were made to adapt the method to various scenarios.First, we tailored the time zone configurations for each site to ensure that DTDS correctly represents the growing and decaying period of the boundary layer according to the geographical location.Second, SGP blind zone is ∼0.2 km, which makes DTDS retrieve the PBLH above this level.We adjusted the blind zone of MAO, TMP, and COR to ensure peak efficiency.Third, instead of using the traditional morning radiosonde measurements, we opted to use the LCL as the upper boundary for determining the initial position of PBLH to avoid the interference of the morning residual layer.The LCL was estimated from surface meteorological data (relative humidity, temperature, pressure) based on an exact expression (Romps, 2017).Therefore, our primary use of the radiosonde data was for validation purposes only.Finally, we expanded the data set's span, incorporating over 20 years of data, marking a significant augmentation from the earlier evaluation duration of 8 years, as reported by Su et al. (2020).
The modifications made here do not alter the original algorithm, as the parameters and architecture remain unchanged.Figure S1 in Supporting Information S1 compares the DTDS-PBLH with and without adjustments to the blind zone.With MAO's lidar blind zone located at lower altitudes than the SGP's blind zone, tweaking the blind zone in DTDS enables more accurate retrieval of PBLH, particularly in shallow boundary layers.It is worth emphasizing that such changes are not to the algorithm (at least not to its physical foundation and values) but to the condition under which it is applied, thus enhancing its potential.

ARM PBLH Product Based on Ceilometer Measurements
ARM provides continuous estimates of the PBLH using ceilometers (i.e., Ceil-PBLH) at four of the five ARM observatories used in this study: SGP, MAO, COR, and TMP.Specifically, ARM deploys Vaisala CL31 ceilometers, which have a maximum vertical range of 7.7 km (Munkel & Rasanen, 2004).The Vaisala CL31 ceilometer provides a total attenuated backscatter coefficient profile at 910 nm, with a vertical resolution of 10 m and a temporal resolution of 2 s (Münkel et al., 2007).The ARM Ceil-PBLH method starts by averaging the backscatter data to a temporal resolution of 16 s, after which the well-known gradient method is applied.This method searches for local gradient minima of the corrected total backscatter signal.Additionally, the ARM Ceil-PBLH product excludes cloud and precipitation cases and suppresses false layer identification (D.Zhang et al., 2022).In this study, we compared the performance of the DTDS algorithm against this existing ARM-PBLH product.

Overview of the DTDS-Derived PBLH Product
Figure 2 shows typical examples of normalized backscatter signals measured by MPLs at the five ARM observatories.Also shown are DTDS-derived PBLHs, radiosonde estimates derived following S. Liu and Liang (2010), and ceilometer-derived PBLHs (based on the gradient method).Figure 2a shows the evolution of the PBL at the SGP in wintertime as an example.During this season, it is common to have little diurnal variation in the boundary layer, so the aerosol structure tends to be highly stratified.This stratification usually complicates the retrieval of the PBLH based on lidar techniques, with most retrieval algorithms performing poorly.In this SGP case, the DTDS algorithm correctly captures the boundary layer's evolution and matches the RS-PBLH estimates, even in the presence of a stratified aerosol layer.The Ceil-PBLH retrievals also correctly captures the PBLH evolution with few jumps between different atmospheric layers.
The MAO example shown in Figure 2b shows a nocturnal boundary layer initially capped with boundary-layer clouds.Clouds induce a strong step signal in the lidar backscatter and depending on the type of cloud, can either facilitate or interfere with the PBLH retrieval.In the MAO example, these boundary-layer clouds help define the PBLH.The DTDS-PBLH and Ceil-PBLH retrievals agree well and are consistent with the RS-PBLH estimates.Around 14 LT, the cloud signal disappears in the backscatter profile, exposing a more complex aerosol structure with multiple aerosol layers and complicating the PBLH retrieval.
While the Ceil-PBLH retrievals jump between different atmospheric layers, sometimes selecting elevated aerosol layers as the PBLH by mistake, the DTDS retrievals consistently track the PBLH.This consistent tracking of the PBLH is due to the selection scheme of DTDS (cf. Figure 3b in Su et al., 2020) and its consideration of the diurnal variability of the boundary layer.DTDS determines whether the boundary layer is decaying, growing, or in other periods and prioritizes the selection of future retrievals closer to current retrievals.In the case of Figure 2b, the strong signal of the PBLH caused by the presence of boundary layer clouds before 14 LT helps DTDS compute the PBLH after 14 LT.The selection scheme would prioritize heights nearer to the cloud's signal in a decaying boundary layer situation.Figure 2c shows a situation similar to MAO during the morning hours, where the PBLH in COR is capped with boundary-layer clouds that accurately define the PBLH.In this example, both the DTDS-PBLH and Ceil-PBLH retrievals are similar throughout the day.Figure 2d shows an example for TMP.In this case, the PBL is also capped with boundary-layer clouds.However, the Ceil-PBLH retrievals considerably fluctuate between different atmospheric levels.One explanation for such a poor performance might be that the ceilometer algorithm computes the PBLH at each time step without considering the previous estimate, likely causing a significant variation of the retrieval at every time point.Note that the DTDS-PBLH algorithm consistently tracks the development of the PBLH without any significant fluctuations in the retrieval.
Unlike the other sites analyzed, the MPL deployed at the FKB site was not fully calibrated, having neglected the afterpulse correction (Figure 2e).This complicates PBLH retrievals from that site, more prominently noticeable under clear-sky conditions.In the example shown, the PBL capped with boundary-layer clouds stays at a relatively constant height (∼1 km) throughout the day.Even without the afterpulse correction, DTDS-PBLH and RS-PBLH retrievals match well.
In addition to examining cases on specific days, it is also insightful to investigate the climatology of the diurnal cycle when assessing the performance of the DTDS algorithm at the five different sites.Figure 3 shows box-andwhisker plots for the diurnal cycles of the DTDS-PBLH and RS-PBLH retrievals at the five ARM sites.In general, the PBLH is a shallow layer during the morning hours between 06:00 to 09:00 a.m.(6-9 LT).The lowest morning layer is found at FKB, where the entire distribution between 6:00 a.m. and 7:00 a.m.(6-7 LT) is below 500 m above ground level (m.a.g.l).COR has the maximum PBLH in the morning hours, with layers extending as high as 1.5 km a.g.l.COR's high topography is one possible reason for this higher PBLH in the morning hours.Since COR is located at 1,141 m, it has a much larger surface roughness and might act as an elevated heat source, both causing strong turbulence and convective instability that favors deep boundary layers.COR's topography also favors the formation of orographic boundary layer clouds, which are frequently observed in this site.Many of these clouds grow to congestus depths.Determining the PBLH for this cloud regime is challenging, and most PBLH retrieval methods would estimate it at the cloud base height.Addressing the challenges in determining the PBLH for cloud regimes that grow to congestus depths requires further and detailed examination, including comprehensive observational data from multiple sources such as Doppler lidars, radiosondes, and cloud radars to accurately identify cloud profiles and the vertical turbulent structure of the PBL.Additionally, large eddy simulations can be used to elucidate the interactions between the PBL and cloud regimes, offering a more complete understanding of the underlying processes.
Independent of the site, the PBLH grows between 8:00 a.m. and 10:00 a.m.(8-10 LT) as the turbulent eddies progressively become larger, responding to a warmer surface.Looking at the mean values (orange lines), this increment is steeper for MAO, a tropical region near the equator with a continuous energy surplus.In this site, the large surface sensible heat flux together with a high surface temperature result in stronger boundary layer turbulence and a relatively faster growing-boundary layer, especially during shallow cumulus days (Tian et al., 2021).FKB also has a steep increment after 8:00 a.m.One of the reasons could be that FKB's measurements were made in the summertime to understand the intense orographic convection in this region.The PBLH reaches a maximum in the afternoon, occurring at different times depending on the site.For instance, the PBLH over the SGP with its extensive and flat terrain is mainly influenced by surface heating, which causes the PBL to grow throughout the day, peaking in the mid to late afternoon (S.Liu & Liang, 2010), consistent with our retrievals at SGP, where the PBLH peaks around 16-17 LT.In contrast, GOAMAZON and CACTI usually experience an earlier peak, as these regions are convectively active (Tian et al., 2021).Frequent convective precipitation in the afternoon can notably suppress the growth of the PBL.The PBLH decays in the late afternoon.
Some of the differences observed in the climatological pattern of each site may also arise from the different data set lengths.The SGP data set spans almost 20 years, allowing this site to have a robust representation of variations related to seasonal forcing, multiple phases of the El Niño-Southern Oscillation ("ENSO"), and other climate oscillations.On the contrary, the measurements made in the Mobile ARM sites correspond to relatively short campaigns with a data set spanning almost 2 years in the best case at MAO.Therefore, the PBLH evolution observed in Figure 3 for the mobile facilities could be highly influenced by unstable and deeper boundary layers occurring in summer or stable and shallower boundary layers in winter.For instance, the FKB site in Germany does not have measurements for the winter months of January, February, and March, potentially leading to an overestimation of the climatology of PBLH for that site.

Evaluation of the DTDS-PBLH Product
Figure 4 compares PBLH estimates obtained by the ceilometer method and derived from radiosondes at four of the five ARM sites.Note that FKB is not included in this figure because the Ceil-PBLH product is not available.Kernel density estimates (KDEs), which is a non-parametric technique that provides a continuous estimate of the probability density function of a data set, in this case, the PBLH, are shown.Also given are the correlation coefficients (Rs), the Root-Mean-Square Errors (RMSEs), and the Mean Absolute Errors (MAEs) between the DTDS-PBLH and the RS-PBLH retrievals.In general, the relation between Ceil-PBLH and RS-PBLH is highly dispersed, with R ranging from 0.48 at SGP to 0.72 at MAO, RMSE oscillating between 0.48 and 0.7 km, and MAE between 0.28 and 0.46 km.
Figure 5 illustrates the same analysis as Figure 4, but instead of using the Ceil-PBLH data set, DTDS-PBLH, and RS-PBLH retrievals are compared.The DTDS-PBLH algorithm outperforms the Ceil-PBLH method at every site, with higher correlations and smaller errors.Figure 4 also shows a strong linear relationship between both sets of retrievals, with R values ranging from 0.77 at COR to 0.93 at MAO.These higher correlations, in comparison with the ceilometer method, attest to the effectiveness of the DTDS algorithm in resolving the temporal evolution of the PBLH.The robustness of the DTDS algorithm is also echoed in the error values.Specifically, RMSEs oscillate between 0.27 and 0.41 km, approximately 0.25 km smaller than the RMSE found with Ceil-PBLH, and MAEs oscillate between 0.18 and 0.27 km.approximately 0.14 km smaller than Ceil-PBLH.Figure 6 summarizes the information from Figures 3 and 4.
Figure 5 shows that the DTDS algorithm tends to overestimate shallow boundary layers, typically below 1 km, and underestimate deep boundary layers.Such a situation is particularly noticeable at the SGP sites.At the MAO site, the overestimation of the PBLH is more noticeable for deep boundary layers.In particular, MAO is located in a tropical region where the PBL expands to significant altitudes during the afternoon, fueled by the surplus of energy prevalent in these regions.However, the DTDS algorithm performs the best at MAO, with the highest R and the smallest error values, surpassing even those of the SGP site, where the DTDS algorithm was developed and tested.
The overestimation of shallow boundary layers generally occurs during the morning hours when the PBL tends to be stable and highly stratified.Both conditions pose difficulties in estimating the PBLH based on backscatter information, which is more related to lidar limitations than to problems in the DTDS algorithm.For example, in the case of extratropical sites like SGP, the performance of the DTDS algorithm varies with season.Retrievals from the SGP site have higher errors during winter because the PBLH tends to be shallower than in other seasons.
The PBLH can sometimes be below the near-surface 150-m MPL blind zone of the region, impeding the correct estimation of the PBLH.

Analysis of Errors and Limitations
The analysis of errors and limitations in the performance of the DTDS algorithm reveals that attributing all inaccuracies solely to the algorithm itself would be overly simplistic.Instead, this subsection comprehensively examines the multiple sources of errors, including those resulting from lidar measurements and limitations, radiosondes measurements, and the RS-PBLH algorithm.By delving into these factors, light will be shed on the intricate interplay between these components, highlighting the inherent complexity of the PBLH retrieval.At the same time, this error analysis will help recognize areas for potential enhancement and underscore the need for further refinement of the DTDS algorithm.

Errors in the Radiosonde Estimation
Despite its limitations, the most reliable and accurate method for measuring PBLH is the use of radiosonde data.The estimations derived from radiosonde PBLH measurements are based on the thermodynamic profile of the atmosphere, which provides valuable information about the stability and structure of the boundary layer.The radiosonde technique is particularly useful in detecting the presence of inversion layers, which play a crucial role in determining the PBLH.In contrast, MPL-based estimations rely on the idea that aerosols are well mixed inside the boundary layer, creating a strong step signal between the PBL and the free atmosphere.The difference  RS-PBLH retrievals are widely accepted in the scientific community as the ground truth.However, even though radiosonde measurements are the most reliable source to estimate the PBLH, the computed PBLH can sometimes be ambiguous.Radiosonde-based thermodynamic profiles can be vague, for example, when there are multiple inversion layers or no inversions.Both situations complicate the PBLH retrieval, leading to great uncertainty.Furthermore, algorithms that use radiosonde information to compute the PBLH are also susceptible to errors, given the multiple parameters needed that are not necessarily adjusted to local environmental conditions, such as the multiple threshold values needed in the Liu and Liang method.The computed PBLH based on radiosonde measurements can also vary depending on the chosen thermodynamic variable, namely the Potential Temperature (θ) or the Virtual Potential Temperature (θ v ) (S. Liu & Liang, 2010).
Figure 7b illustrates the PBLH computed at MAO on 14 September 2015, along with the backscatter information during a relatively clear-sky day with a highly stratified aerosol profile.Note the elevated aerosol layer signal during the early morning that progressively disappears as the day evolves.At noon, the DTDS-PBLH algorithm retrieves the PBLH around 1.5 km, corresponding to the approximate location of the aerosol layer.The RS-PBLH algorithm identifies the PBLH at ∼2.5 km.To better understand this discrepancy, Figure 7a presents θ and θ v profiles measured near noon.Also shown is the height of the PBLH identified by the RS-PBLH algorithm.In both temperature profiles, there is one capping inversion near 1.5 km and a weaker one around 2.5 km, where the Liu and Liang algorithm estimates the PBLH.Multiple inversion layers in a temperature profile pose a challenge when attempting to identify the PBLH based solely on thermodynamic considerations because deciding which of the inversion layers should be considered as the PBLH can be ambiguous.In this case, lidar data provides extra information to better describe the PBL.Here, the capping inversion around 1.5 km aligns with the DTDS-PBLH estimate and the extension of the aerosol layer, suggesting that this position would be more appropriate as the PBLH instead of the second inversion around 2.5 km, which the RS-PBLH algorithm assigns as the PBLH.
While the Liu and Liang method is generally reliable in estimating the PBLH under a variety of thermodynamic conditions, it encounters a challenge in identifying the inversion layer in this specific scenario due to the predetermined threshold values (δ s , δ u , and θr ).These values, derived empirically in S. Liu and Liang (2010), may lead to erroneous estimations in certain situations, such as when multiple inversion heights are present.However, it's important to note that the method's overall solid performance remains intact.In Figure S2 in Supporting Information S1, we present a thorough sensitivity analysis of the effect of the threshold values on PBLH determination based on the MAO data set.This analysis reveals that varying threshold values may result in shallower or deeper boundary layers compared to the standard values, highlighting the importance of these threshold values on the overall algorithm's performance.Despite the variations observed, the Liu and Liang method consistently provides reliable PBLH retrievals, and using the standard values helps with reproducibility and comparison purposes.
As mentioned previously, differences between the RS-PBLH and the DTDS-PBLH retrievals also depend on the type of temperature used to compute the PBLH. Figure 9 is like Figure 5 but uses θ v instead of θ to compute the RS-PBLH.While these figures are similar, they are not identical.Overall, the PBLH computed using θ performs better than the one using θ v because θ v contains considerable uncertainties due to errors measuring humidity (S.Liu & Liang, 2010).However, for three of the five sites, the correlation coefficient between RS-PBLH and DTDS-PBLH retrievals is higher, going from 0.77 to 0.79 at SGP, from 0.85 to 0.87 at TMP, and from 0.78 to 0.81 at FKB, albeit with a slightly higher RMSE of ∼0.1 km.
The higher correlation observed at SGP and TMP when using θ v is attributed to its superior representation of boundary layer thermodynamics.θ v accounts for the effects of moisture content in the air, whereas θ does not.Since water vapor is less dense than dry air, moist, unsaturated air is more buoyant than dry air of the same temperature.Therefore, θ v better represents buoyancy-driven motions in the atmosphere than θ does and, in principle, should be the preferred variable for PBLH estimation if errors in humidity measurements were small, but in reality this is not the case.These differences in the RS-PBLH retrievals arising from the temperature used reveal that the dispersion observed in the KDE functions of Figures 5 and 9 are influenced not only by the performance of the DTDS algorithm but also by the various complexities inherent in the RS-PBLH.This finding suggests that there is room for improvement in the accuracy of RS-PBLH estimates as well.
The distinction between RS-PBLH calculated using θ and θ v is demonstrated in Figure 10, showing KDEs between RS-PBLH based on θ and θ v at the five ARM sites.Among these sites, the two sets of RS-PBLH estimates agree the best at FKB, with an R of 0.95 and an RMSE of 0.2 km.The worse agreement is seen at MAO, with an R of 0.92 and an RMSE of 0.33 km.At all sites except FKB, the RS-PBLH based on θ v consistently yields higher values than the RS-PBLH based on θ, contributing significantly to the larger RMSE observed, particularly at locations like MAO.
Finally, it's important to emphasize that our discussion on the limitations and potential sources of errors in the Liu and Liang technique is not intended to cast any doubt on the overall quality of radiosonde-derived PBLH, which remains the standard method.Rather, we aim to explore potential mismatches between the MPL and radiosonde techniques and highlight areas for further investigation and refinement in atmospheric boundary layer studies.

Errors and Limitations of the MPL
Another important source of errors arising when computing the DTDS-PBLH is that resulting from lidar limitations.For example, Figure 11 showcases a situation that occurred over TMP on 15 February 2014.On that day, the presence of a stratocumulus-topped boundary layer (STBL) attenuated the MPL backscatter, preventing the instrument from sensing the atmosphere beyond the cloud-base height.In this case, the DTDS algorithm retrieves the PBLH at the same height as the cloud base at ∼0.7 km.However, when looking at the temperature profiles around noon in Figure 11a, a well-mixed boundary layer with a strong capping inversion at ∼1.1 km is seen.This inversion aligns with the cloud-top height shown in Figure 11a.Considering that the PBLH is defined as the maximum altitude influenced by the earth's surface, the RS-PBLH would be a better estimate of the PBLH.The presence of an STBL and the limited capability of the MPL under cloudy conditions introduces errors in such cases.
The performance of the DTDS algorithm is notably diminished in instances where the lidar instrument is not operating properly.The COR site provides an example of the errors that arise from poor lidar performance.The initial MPL installed at the COR site was misaligned, leading to a consistent spurious signal always present in the backscatter profile.Figures 12a and 12b shows two different days when the MPL at COR was misaligned, one under clear-sky conditions and a second under cloudy conditions.In both cases, a strong constant step signal is visible between 0.5 and 1 km.For reference, this spurious signal is not observed in the backscatter signal of the ceilometer (see Figures 12c and 12b).The spurious signal in the backscatter profile of the MPL impeded the accurate determination of the PBLH on the clear-sky day, resulting in the DTDS algorithm incorrectly placing the   PBLH at ∼0.8 km.In the case when boundary-layer clouds were present (Figure 12b), the strong backscatter signal from the clouds counteracted the spurious backscatter signal, helping the DTDS algorithm accurately retrieve the PBLH and resulting in a better agreement between the DTDS-PBLH, Ceil-PBLH, and RS-PBLH retrievals.
To address the misalignment issue in the MPL, a new lidar was installed at the site on 18 February 2019.This replacement successfully eliminated the spurious signal, leading to a notable improvement in the performance of the DTDS algorithm.Figure S3 in Supporting Information S1 shows an example of the MPL backscatter signal after the sensor was replaced.Figure 13 compares the DTDS performance at retrieving the PBLH before and after the MPL was changed under clear-sky and cloudy conditions.Note that erroneous DTDS-PBLH estimates were included because the primary aim is to highlight the pronounced improvements observed in the DTDS performance exclusively resulting from replacing the lidar.Overall, R is higher (Figure 13a), and the RMSE is smaller (Figure 13b) under cloudy conditions than under clear-sky conditions.Also, after the MPL was replaced, R increased from 0.58 to 0.7 on average, and RMSE decreased from 0.48 to 0.37 km.In the case of the MAEs, we also see a slight improvement when the lidar is changed, even considering that the MAE weighs all errors equally.

Conclusions
Gaining a comprehensive understanding of the PBLH is crucial for advancing our knowledge across numerous atmospheric processes and applicable fields.Despite the longstanding reliance on radiosondes for PBLH measurements, there is a pressing need to enhance the precision of PBLH estimates, particularly when applying lidar techniques.Over recent years, lidar systems have been increasingly used to estimate the PBLH, primarily because they overcome a key limitation of traditional radiosonde estimates, namely, the lack of temporal resolution.Lidar methods allow for continuous PBLH monitoring, thereby advancing our grasp of atmosphere-surface processes and the evolution of the PBL.However, lidar systems have their own inherent limitations to incur considerable uncertainties and biases, especially under stable and highly stratified conditions.
The DOE ARM project has deployed an extensive network of MPL systems in many places, thereby presenting a valuable opportunity to investigate PBLH in different parts of the world.In this study, we took advantage of this large network by evaluating the performance of the DTDS algorithm developed by our team (Su et al., 2020).It was aimed to tackle with a common problem of incoherent variations between atmospheric thermodynamics and the backscattering of aerosols chiefly in early morning by modifying the traditional wavelet method for enhancing temporal and spatial coherence.The DTDS algorithm was developed and tested for calculating PBLHs at the SGP site, the ARM's primary permanent sites where it shows superior performance over the conventional wavelet and gradient methods attested by higher correlations and smaller RMSE but its performance at other places of different conditions are unknown.
The DTDS is thus evaluated at five ARM observatories situated in different continents (North America, South America, and Europe) where more than 6 months of MPL data were available.Notably, the DTDS-PBLH and the radiosonde-derived PBLH agreed well, with R ranging from 0.77 at the COR site to 0.93 at the MAO site, and the root-mean-square error ranging 0.27 km at MAO to 0.41 km at SGP.The correlation and the RMSE instill confidence in the new method's robustness in representing the boundary layer's temporal variability, thus unveiling new possibilities for studying surface atmospheric processes.Compared to existing PBLH products for the same sites, for example, Ceil-PBLH, the RMSEs of DTDS-PBLH were roughly a third smaller than the RMSE of Ceil-PBLH.On the other hand, the MAEs of DTDS-PBLH were half smaller than the MAEs of Ceil-PBL.Current lidar-based RMSE of PBLH estimates at the five different sites were close to 0.6 km, which is reduced to less than 0.4 km for the DTDS-derived PBLH.In the case of the MAE, the current ceil-PBLH errors were close to 0.41 km that is reduced to an average of 0.23 km for the DTDS.
Despite the relatively superior performance, both the DTDS algorithm and input data suffer from limitations, leading to retrieval errors in the PBL product.Besides, there also exist uncertainties in radiosonde-based PBLH estimates due to the complexity of atmospheric thermodynamics.For example, at the COR site, lidar misalignments in the first MPL installed there produced a spurious signal that hampered the DTDS performance.Further limitations in the DTDS algorithm include constraints in the MPL visibility, such as during cloudy conditions that attenuate the backscatter signal.Radiosonde measurements may also introduce ambiguity in determining the PBLH, causing disparities between radiosonde and DTDS estimates, especially when multiple inversion layers are present or when there are no inversions.Moreover, the thermodynamic variable used in calculating the radiosonde-derived PBLH (e.g., potential temperature or virtual potential temperature) can introduce variability in the results.A thorough understanding of these errors and limitations is pivotal for continuous improvements in the DTDS algorithm, enhancing its performance and reliability in estimating PBLH.

Figure 1 .
Figure 1.Geographical distribution of the five Department of Energy Atmospheric Radiation Measurement sites used in this investigation.

Figure 2 .
Figure 2. Examples of the computed planetary boundary layer height at the five Atmospheric Radiation Measurement (ARM) sites using the Different Thermo-Dynamic Stabilities algorithm (black dots), the ARM ceilometer product (white dots), and radiosonde estimates (magenta stars).The colormap corresponds to the normalized backscatter signal measured by the Micropulse Lidar at each site.

Figure 3 .
Figure 3. Climatologies of the diurnal cycle of planetary boundary layer height (PBLH) computed using the Different Thermo-Dynamic Stabilities algorithm-PBLH and radiosondes (RS-PBLH) at the five Atmospheric Radiation Measurement observatories.The horizontal orange lines, blue bars, and whiskers represent median values, interquartile ranges, and ranges of the data, respectively.Time along the x-axis is Local Time.

Figure 4 .
Figure 4. Kernel distribution estimates of Atmospheric Radiation Measurement (ARM) ceilometer-and radiosonde-derived planetary boundary layer heights at the five ARM sites.The correlation coefficient (R), root-mean-square error, and mean absolute error are given in each panel.Dashed red lines correspond to 1:1 lines.

Figure 5 .
Figure 5. Kernel distribution estimates for the micropulse lidar-based Different Thermo-Dynamic Stabilities-planetary boundary layer height (DTDS-PBLH) and the radiosonde PBLH based on the Liu-Liang method using virtual potential temperatures at the five Atmospheric Radiation Measurement sites.The correlation coefficient (R), root-mean-square error, and mean absolute error are given in each panel.Dashed red lines correspond to 1:1 lines.

Figure 6 .
Figure 6.Correlation coefficients, root-mean-square errors, and mean absolute errors from the Different Thermo-Dynamic Stabilities-and radiosonde-derived planetary boundary layer height (PBLH) relations (red bars) and the ceilometer-and radiosonde-derived PBLH relations (blue bars) at the five atmospheric radiation measurement sites.

Figure 7 .
Figure 7. Planetary boundary layer heights (PBLHs) at MAO on 14 September 2015.(a) Profiles of potential temperature (θ) and virtual potential temperature (θ v ) measured around noon on the same date (black and blue lines, respectively).The magenta line corresponds to the RS-PBLH.In this example, RS-PBLHs based on θ and θ v , and the Atmospheric Radiation Measurement RS-PBLH have the same values.Therefore, only one line is shown.(b) Backscatter information for the same day where the colormap indicates the normalized backscatter signal.Black dots correspond to Different Thermo-Dynamic Stabilities-PBLHs, while stars are RS-PBLHs.
Figure 8 provides a sensitivity analysis of the threshold values for the same date shown in Figure 7.The red lines represent the PBLH using the standard values, while the blue line indicates a modified version.As demonstrated by adjusting the standard values of θr and δ u , the Liu and Liang algorithm identifies the first inversion layer at approximately 1.5 km as the PBLH, aligning with the PBLH-DTDS retrieval.

Figure 8 .
Figure 8. Sensitivity analysis comparing the threshold values proposed by Liu and Liang for (a) Θ r , (b) δ s , and (c) δ u at MAO on 14 September 2015.The red lines represent the standard threshold values, while the blue line depicts the modified thresholds.

Figure 9 .
Figure 9. Same as Figure 5, but instead of using θ to compute the planetary boundary layer height, θ v was used.

Figure 11 .
Figure 11.Planetary boundary layer heights (PBLHs) at TMP on 22 October 2014.(a) Profiles of potential temperature (θ) and virtual potential temperature (θ v ) measured around noon on the same date (black and blue lines, respectively).The magenta line corresponds to the RS-PBLH.For this example, RS-PBLH based on θ and θ v , and the Atmospheric Radiation Measurement RS-PBLH have the same values.Therefore, only one line is shown.(b) Backscatter information for the same day where the colormap indicates the normalized backscatter signal.Black dots correspond to Different Thermo-Dynamic Stabilities (DTDS-PBLH), stars are RS-PBLHs, white dots are cloud-top heights (CTH), and gray dots are cloud-base heights (CBH).

Figure 10 .
Figure 10.Kernel distribution estimates for radiosonde planetary boundary layer heights based on potential temperature and virtual potential temperature at the five Atmospheric Radiation Measurement sites.The correlation coefficient (R), rootmean-square error (RMSE), and mean absolute error (MAE) are given in each panel.Dashed red lines correspond to 1:1 lines.

Figure 12 .
Figure 12.Planetary boundary layer heights (PBLHs) at COR using the Different Thermo-Dynamic Stabilities algorithm (black dots), the Atmospheric Radiation Measurement Ceil-PBLH product (white dots), and radiosonde estimates (magenta stars) prior to the replacement of the micropulse lidar (MPL) on (a, c) a clear-sky day and (b, d) a cloudy day.The colormap indicates the normalized backscatter signal measured by the (a, b) MPL and the (c, b) ceilometer.The constant orange-red belt observed in (a, b) between 0.5 and 1 km represents spurious signals.

Figure 13 .
Figure 13.(a) Correlation coefficients, (b) root-mean-square errors, and (c) mean absolute errors from the DTDS-and radiosonde-derived planetary boundary layer height relations before and after the micropulse lidar was replaced at the COR site.Red bars correspond to clear-sky conditions, and blue bars represent cloudy conditions.