A Comparison of Wintertime Atmospheric Boundary Layer Heights Determined by Tethered Balloon Soundings and Lidar at the Site of SACOL

: High-precision and -resolution atmospheric boundary layer height (BLH) has received increasing attention in air pollution research in recent years. The low time resolution of sounding data is the main challenge to validate BLH retrieval from lidar observations. To resolve this issue, we conducted simultaneous tethered balloon sounding and lidar observations at the Semi-Arid Climate and Environment Observatory of Lanzhou University (SACOL) during winter 2019–2020. The BLHs derived from the tethered balloon sounding data were 170, 210, 393, 676, 423, and 190 m at 02:00, 08:00, 11:00, 14:00, 17:00, and 20:00 (Beijing time), respectively. The diurnal evolution of BLH was reasonably captured by lidar observation-based wavelet covariance transform and ideal proﬁle ﬁtting methods, which exhibited correlation coefﬁcients of 0.91 and 0.89, respectively, with the BLHs determined from tethered balloon sounding data. The lidar results slightly overestimated the BLHs, though all results were acceptable when considering both the absolute and relative errors with respect to BLHs from the tethered balloon data. Our results revealed high-precision and -resolution diurnal variations in BLH at SACOL in Northwest China and suggest the importance of validating lidar-based BLHs using simultaneous sounding data. Our results revealed high-resolution and -accuracy diurnal variations of BLHs based on tethered balloon sounding observations at the site of SACOL in Northwest China. The diurnal variation in retrieval errors in lidar-BLHs suggests the importance of validating lidar-based BLHs using simultaneous sounding data. To our best knowledge, our study presents the ﬁrst comprehensive high-resolution results (3 h time resolution during daytime) on the comparison between lidar-BLHs and TBS-BLHs in China.


Introduction
The atmospheric boundary layer is the lowest layer in the atmosphere, directly forced by the ground, with obvious turbulence characteristics [1]. Energy and material exchange between the earth and atmosphere occurs in this layer. Thermal and dynamic forcing from the ground causes characteristic diurnal variations of the atmospheric boundary layer structures, which directly affects the diffusion and transportation of air pollutants [2][3][4] and indirectly impacts human health and productivity [5]. The atmospheric boundary layer height (BLH), generally defined as the height where turbulence almost disappears, determines the atmospheric environmental capacity [6]. The concentration of particulate matter <2.5 µm in diameter (PM 2.5 ) is reversely correlated with BLH [7][8][9][10], and BLH contributes 25% of total PM 2.5 in Eastern China [11].
A mobile depolarization lidar (MiniMPL; Sigma Space, USA; Figure 1e) was installed 50 m away from the tethered balloon to observe atmospheric aerosols at 532 nm with a blind zone of 100 m, a vertical distribution of 30 m, and a time resolution of 1 min (Table S2). The calculation of normalized range-corrected backscattering signal (NRB) was as follows: where Raw(z) is the raw lidar data (counts µs −1 ), D(z) is the dead time correction, AP(z) is the afterpulse correction (Hz), B is the background signal correction (counts km 2 µs −1 ), O(z) is overlap correction, and E is the laser energy (µJ). A total of 18 days of lidar observations without clouds and precipitation were used to retrieve BLH.

Determining BLH from TBS Data
The structure of atmospheric boundary layer has diurnal variation, in which convective boundary layers predominate during daytime and stable boundary layers at night, and it is also complex during the transition between daytime and nighttime. The convective boundary layer height was determined from TBS data in the daytime as: (1) the height at which virtual potential temperature (θ v ) is the same as the surface value (Holzworth method or parcel method) [43,59] and (2) the height at which the bulk Richardson number (Rb) is the same as 0.25 (named Rb number method) [44,60]. The bulk Richardson number was calculated as: where θ v is the virtual potential temperature, g is the gravitational acceleration, u and v are wind speed components, b is a constant, and u * is the surface friction velocity, which are so much smaller than the bulk shear terms in the denominator that they can be ignored [16,17]. The z is the height, and s represents the surface (2 m). Here, u s and v s were set to u and v at height of 2 m. This setting is further investigated in the Supplementary Materials ( Figure S1). The correlation coefficients between BLHs derived from Holzworth and Rb number method were 0.99 for the matched 26 daytime profiles ( Figure S2a). The Holzworth method can be understood as a simplification of the Rb number method where the shear contribution is neglected. Thus, they are only applicable for unstable conditions [44]. It is difficult to derive nocturnal stable boundary layer height because of the influence of gravity waves, intermittent turbulence, and other complex motions [44]. The stable boundary layers can be derived from temperature and wind profiles [61][62][63]. However, more than half of the wind shear method-derived BLHs were lower than 100 m while most of the temperature inversion-derived BLHs varied from 100 m to 300 m ( Figure S3). Seidel et al. (2010) indicated that a surface-based temperature inversion is a clear indicator of the stable boundary layer [45]. Therefore, stable boundary layer height was determined from TBS data in the nighttime as the height at the top of a surface-based temperature inversion. For comparison, the Rb number method was also used to derive stable boundary layer height [64,65]. Here, us and vs were set to zero to avoid extremely low BLHs ( Figure S4). The correlation coefficients between BLHs derived from temperature inversion and Rb number method were 0.60 for the matched 30 nighttime profiles ( Figure S2b).
After comparison, the Holzworth method and temperature inversion method were selected to derive the BLHs based on TBS (TBS-BLHs) in daytime and nighttime, respectively ( Figure 2). Moreover, the finer structures of the atmospheric boundary layer were analyzed by manual distinction with some principles (in Section 3.1): (1) the residual capping inversion layers are the thin layers beneath the free atmosphere with a temperature inversion; (2) the residual mixed layers are above the nocturnal stable boundary layers and below the residual capping inversion layers; (3) the daytime capping inversion layers (also referred to as entrainment zones) are the thin layers above the convective boundary layers with a temperature inversion. capping inversion layers are the thin layers beneath the free atmosphere with a temperature inversion; (2) the residual mixed layers are above the nocturnal stable boundary layers and below the residual capping inversion layers; (3) the daytime capping inversion layers (also referred to as entrainment zones) are the thin layers above the convective boundary layers with a temperature inversion.

Determining BLH from Lidar Observations
We applied the WCT [34,66] and FIT [35,67] methods to determine BLHs from lidar observations as these methods have been widely used in previous studies [42,55]. Since atmospheric aerosols are generally retained in the atmospheric boundary layer, the backscatter signal always decreases rapidly in the transition zone from the atmospheric boundary layer to the free troposphere and the WCT can amplify the gradient in the transition zone [68], making it easier to determine the BLH. Specifically, the wavelet basis function was defined as: where b is the wavelet translation, a is the wavelet spatial dilation, and z is the backscatter signal profile height. The wavelet covariance transform function wf(a, b) was defined as: where b(z) is the backscatter signal at height z and zb and zt are the lower and upper height limit. At a given value of a, the BLH is defined as the value of b corresponding to maximum wavelet covariance when b steps from zb to zt (Figure 3a). The value of a impacts the BLHs of the WCT method to some extent [41]. The value of zb should be above the height of a/2 (Equation 2), which requires a small a value for shallow atmospheric boundary layer under conditions such as nocturnal stable atmospheric boundary layer in winter. On the other hand, a relatively large a value should be

Determining BLH from Lidar Observations
We applied the WCT [34,66] and FIT [35,67] methods to determine BLHs from lidar observations as these methods have been widely used in previous studies [42,55]. Since atmospheric aerosols are generally retained in the atmospheric boundary layer, the backscatter signal always decreases rapidly in the transition zone from the atmospheric boundary layer to the free troposphere and the WCT can amplify the gradient in the transition zone [68], making it easier to determine the BLH. Specifically, the wavelet basis function was defined as: where b is the wavelet translation, a is the wavelet spatial dilation, and z is the backscatter signal profile height. The wavelet covariance transform function wf(a, b) was defined as: where b(z) is the backscatter signal at height z and z b and z t are the lower and upper height limit. At a given value of a, the BLH is defined as the value of b corresponding to maximum wavelet covariance when b steps from z b to z t (Figure 3a). The value of a impacts the BLHs of the WCT method to some extent [41]. The value of z b should be above the height of a/2 (Equation (2)), which requires a small a value for shallow atmospheric boundary layer under conditions such as nocturnal stable atmospheric boundary layer in winter. On the other hand, a relatively large a value should be applied for deep atmospheric boundary layer to minimize noise interference. Hence, in consideration of the diurnal variations of atmospheric boundary layer, the value of a was set to 60 m for samples at 02:00, 180 m at 08:00 and 20:00, and 300 m at 11:00, 14:00, and 17:00, and Compton et al. (2013) also had similar settings [50]. applied for deep atmospheric boundary layer to minimize noise interference. Hence, in consideration of the diurnal variations of atmospheric boundary layer, the value of a was set to 60 m for samples at 02:00, 180 m at 08:00 and 20:00, and 300 m at 11:00, 14:00, and 17:00, and Compton et al. (2013) also had similar settings [50]. FIT defines an ideal profile B(z) as: where Bm is the mean backscatter of the mixed layer, Bu is the mean backscatter in the air immediately above the entrainment zone, zm is a trial value for BLH, and s is related to the depth of the entrainment zone. Based on a simulated annealing algorithm [69], B(z) was fitted to the observed profile b(z) by minimizing root-mean-square deviation between them to determine the four idealized backscatter profile parameters Bm, Bu, zm, and s. Thus, B(z) represents the overall information content of the b(z) and FIT, if free from local extremes within the profile (Figure 3b). The logarithm of the NRB (ln (NRB)) instead of the NRB was used to derive BLH from lidar observations in the present study, which can reduce the underestimation under convective conditions ( Figure S5) [34] and may increase the overestimation in the nighttime. Therefore, we used the ln (NRB) signal in the daytime and NRB signal in the nighttime based on 5 min temporal averages and Perona-Malik anisotropic diffusion filters to improve the signal-to-noise ratio [70]. Finally, the continuity of the BLHs was checked to avoid the influence of transported aerosol layers.

Diurnal Variations in Atmospheric Boundary Layer from TBS Data
Firstly, diurnal variations in atmospheric boundary layer were analyzed based on the TBS data ( Figure 4). The TBS data on 12 December and 24 December 2019 were selected as the typical examples because (1) all the six-time tethered balloon soundings were successfully launched on these two days; and (2) clear sky weather on these two days made it representative for typical atmospheric boundary layer characteristics in winter at SACOL. The surface-based temperature inversion layer was evident in the temperature and virtual potential temperature profiles at 02:00, which was defined as the FIT defines an ideal profile B(z) as: where B m is the mean backscatter of the mixed layer, B u is the mean backscatter in the air immediately above the entrainment zone, z m is a trial value for BLH, and s is related to the depth of the entrainment zone. Based on a simulated annealing algorithm [69], B(z) was fitted to the observed profile b(z) by minimizing root-mean-square deviation between them to determine the four idealized backscatter profile parameters B m , B u , z m , and s. Thus, B(z) represents the overall information content of the b(z) and FIT, if free from local extremes within the profile (Figure 3b). The logarithm of the NRB (ln (NRB)) instead of the NRB was used to derive BLH from lidar observations in the present study, which can reduce the underestimation under convective conditions ( Figure S5) [34] and may increase the overestimation in the nighttime. Therefore, we used the ln (NRB) signal in the daytime and NRB signal in the nighttime based on 5 min temporal averages and Perona-Malik anisotropic diffusion filters to improve the signal-to-noise ratio [70]. Finally, the continuity of the BLHs was checked to avoid the influence of transported aerosol layers.

Diurnal Variations in Atmospheric Boundary Layer from TBS Data
Firstly, diurnal variations in atmospheric boundary layer were analyzed based on the TBS data ( Figure 4). The TBS data on 12 December and 24 December 2019 were selected as the typical examples because (1) all the six-time tethered balloon soundings were successfully launched on these two days; and (2) clear sky weather on these two days made it representative for typical atmospheric boundary layer characteristics in winter at SACOL. The surface-based temperature inversion layer was evident in the temperature and virtual potential temperature profiles at 02:00, which was defined as the nocturnal stable boundary layer. The stable boundary layers occurred at 110 m on December 12 and 160 m on December 24, respectively, at 02:00. The residual mixed layers were defined by tops at 220 and 620 m, respectively. Shallow residual capping inversion layers occurred above the residual mixed layers at 02:00.
nocturnal stable boundary layer. The stable boundary layers occurred at 110 m on December 12 and 160 m on December 24, respectively, at 02:00. The residual mixed layers were defined by tops at 220 and 620 m, respectively. Shallow residual capping inversion layers occurred above the residual mixed layers at 02:00. The stable boundary layer gradually developed upwards with nighttime radiative cooling until sunrise, with tops reaching 170 m on December 12 and 190 m on December 24, respectively, at 08:00. Residual mixed layers developed more rapidly than stable boundary layers from 02:00 to 08:00. After sunrise, surface-based stable boundary layers were destroyed by solar radiative heating and convective boundary layers began to develop, leading to shallow convective boundary layers at 08:00. As solar radiation continued to intensify, convective boundary layers and capping inversion layers continued to develop, with the latter's tops reaching 365 m and 242 m, respectively, at 11:00.
Convective boundary layers reached maximal heights of 910 m on December 12 and 1050 m on December 24, respectively, at 14:00 when turbulent mixing was strong. Then, they weakened, with tops decreasing to 528 and 266 m, respectively, at 17:00. Generally, The stable boundary layer gradually developed upwards with nighttime radiative cooling until sunrise, with tops reaching 170 m on 12 December and 190 m on 24 December, respectively, at 08:00. Residual mixed layers developed more rapidly than stable boundary layers from 02:00 to 08:00. After sunrise, surface-based stable boundary layers were destroyed by solar radiative heating and convective boundary layers began to develop, leading to shallow convective boundary layers at 08:00. As solar radiation continued to intensify, convective boundary layers and capping inversion layers continued to develop, with the latter's tops reaching 365 m and 242 m, respectively, at 11:00.
Convective boundary layers reached maximal heights of 910 m on 12 December and 1050 m on 24 December, respectively, at 14:00 when turbulent mixing was strong. Then, they weakened, with tops decreasing to 528 and 266 m, respectively, at 17:00. Generally, convective boundary layers disappear gradually after sunset and stable boundary layers begin to develop. On 12 December at 20:00, the convective boundary layer disappeared, and a stable boundary layer developed at 20:00, but on 24 December at 20:00 the convective boundary layer was maintained (though very weak). Atmospheric boundary layer evolution was properly characterized by the TBS data. The TBS-BLHs on 12 December The virtual potential temperature profiles at 02:00, 08:00, 11:00, 14:00, 17:00, and 20:00 were averaged to represent wintertime diurnal variations in the atmospheric boundary layer at SACOL (Figure 5), which clearly showed the development of the convective boundary layer and nocturnal stable boundary layer. The average BLHs derived from these profiles were 170, 210, 393, 676, 423, and 190 m at 02:00, 08:00, 11:00, 14:00, 17:00, and 20:00, respectively. The stable boundary layers at 20:00 might be inaccurate because convective boundary layers also contributed to BLHs at this time. The virtual potential temperature profiles at 02:00, 08:00, 11:00, 14:00, 17:00, and 20:00 were averaged to represent wintertime diurnal variations in the atmospheric boundary layer at SACOL (Figure 5), which clearly showed the development of the convective boundary layer and nocturnal stable boundary layer. The average BLHs derived from these profiles were 170, 210, 393, 676, 423, and 190 m at 02:00, 08:00, 11:00, 14:00, 17:00, and 20:00, respectively. The stable boundary layers at 20:00 might be inaccurate because convective boundary layers also contributed to BLHs at this time. More detailed diurnal variations in atmospheric boundary layers can be obtained by using sounding data with high time resolution. Our atmospheric boundary layer diurnal variation results had the highest time resolution published for Northwest China, other than those in Zhang and Wang (2008) [71] based on the eight profiles taken per day. The observation in Zhang and Wang (2008) [71] was conducted twenty years ago during May to June 2000 and the site of the observation was in an arid region of Dunhuang. The influence of atmospheric boundary layer on the surface air pollution concentration in Northwest China was studied using the operational radiosonde sounding data at 08:00 and 20:00 in recent years [72]. The operational radiosonde sounding data were only available during the local early morning (08:00) and late afternoon (20:00) transition times, making it impossible to study the issue under typical convective boundary layer and stable boundary layer conditions in the region. The high-precision and -resolution BLH results in the present study provide valuable information for urban air pollution studies and predications in the study region.

Diurnal Variations in Atmospheric Boundary Layer from Lidar Data
The lidar observations on 12 December and 24 December 2019 were used to assess the diurnal variations of BLH using WCT and FIT methods ( Figure 6). Near the ground, More detailed diurnal variations in atmospheric boundary layers can be obtained by using sounding data with high time resolution. Our atmospheric boundary layer diurnal variation results had the highest time resolution published for Northwest China, other than those in Zhang and Wang (2008) [71] based on the eight profiles taken per day. The observation in Zhang and Wang (2008) [71] was conducted twenty years ago during May to June 2000 and the site of the observation was in an arid region of Dunhuang. The influence of atmospheric boundary layer on the surface air pollution concentration in Northwest China was studied using the operational radiosonde sounding data at 08:00 and 20:00 in recent years [72]. The operational radiosonde sounding data were only available during the local early morning (08:00) and late afternoon (20:00) transition times, making it impossible to study the issue under typical convective boundary layer and stable boundary layer conditions in the region. The high-precision and -resolution BLH results in the present study provide valuable information for urban air pollution studies and predications in the study region.

Diurnal Variations in Atmospheric Boundary Layer from Lidar Data
The lidar observations on 12 December and 24 December 2019 were used to assess the diurnal variations of BLH using WCT and FIT methods ( Figure 6). Near the ground, there were relatively high ln (NRB) signal values during nighttime and low values during daytime, indicating the influence of BLH on the vertical distribution of aerosols. The contribution of aerosol humidification on the lidar signal during the observation period was small due to low relative humidity (42% during nighttime and 51% during daytime) [73,74]. Overall, the vertical distribution of atmospheric aerosols in the atmospheric boundary layer Remote Sens. 2021, 13, 1781 9 of 16 was controlled by the BLH to a large extent and the lidar profiles could determine BLH reasonably. Both WCT and FIT determined BLH locations at transitional places within the lidar profiles and properly captured the diurnal evolution of BLHs. The WCT-BLHs and FIT-BLHs generally matched the TBS-BLHs, though they seemed to be slightly overestimated at 02:00, 08:00, and 20:00 for the selected two days, which may have been caused by aerosols in the elevated residual mixed layers. Thus, nighttime lidar observations should be carefully treated to obtain reasonable BLHs. The WCT and FIT methods generally obtained consistent diurnal variations in BLH, though the latter exhibited smaller fluctuations than the former. there were relatively high ln (NRB) signal values during nighttime and low values during daytime, indicating the influence of BLH on the vertical distribution of aerosols. The contribution of aerosol humidification on the lidar signal during the observation period was small due to low relative humidity (42% during nighttime and 51% during daytime) [73,74]. Overall, the vertical distribution of atmospheric aerosols in the atmospheric boundary layer was controlled by the BLH to a large extent and the lidar profiles could determine BLH reasonably. Both WCT and FIT determined BLH locations at transitional places within the lidar profiles and properly captured the diurnal evolution of BLHs. The WCT-BLHs and FIT-BLHs generally matched the TBS-BLHs, though they seemed to be slightly overestimated at 02:00, 08:00, and 20:00 for the selected two days, which may have been caused by aerosols in the elevated residual mixed layers. Thus, nighttime lidar observations should be carefully treated to obtain reasonable BLHs. The WCT and FIT methods generally obtained consistent diurnal variations in BLH, though the latter exhibited smaller fluctuations than the former. The WCT-BLH and FIT-BLH were comparable (~750 m) but both were lower than the TBS-BLH (1020 m) at 14:00 on December 24 ( Figure 6b). As seen from the TBS temperature profile at 14:00 on December 24 (Figure 4b), a few shallow inversion layers were present in the atmospheric boundary layer, which inhibited the vertical mixing of aerosols. The inversion layers in temperature profile from 500 m to 600 m led to the highest The WCT-BLH and FIT-BLH were comparable (~750 m) but both were lower than the TBS-BLH (1020 m) at 14:00 on 24 December (Figure 6b). As seen from the TBS temperature profile at 14:00 on 24 December (Figure 4b), a few shallow inversion layers were present in the atmospheric boundary layer, which inhibited the vertical mixing of aerosols. The inversion layers in temperature profile from 500 m to 600 m led to the highest gradient in ln (NRB) signal at~750 m, which was defined as the top of the atmospheric boundary layer by both WCT and FIT methods. Moreover, there was a cold front from the north and the wind speed was high on 24 December (Figures S6 and S7) which makes it more complicated to retrieve BLH from both TBS and lidar data. The example on 24 December suggested that a large difference may be exhibited in BLHs from lidar and TBS data; improving the retrieval method of lidar-BLH does not help to yield accurate BLH from lidar observations under such complicated conditions. Both WCT and FIT methods reasonably captured the diurnal evolution of BLH during the observation period at SACOL (Figure 7). The FIT-BLHs exhibited good consistency with the WCT-BLHs but seemed to be slightly higher than WCT-BLHs, especially during the daytime. This overestimation might be attributed to the influence of elevated aerosols in the residual mixed layers. The standard errors of the WCT-BLHs were generally higher than those of the FIT-BLHs, indicating better stability of the FIT method in determining BLH from lidar observations. gradient in ln (NRB) signal at ~750 m, which was defined as the top of the atmospheric boundary layer by both WCT and FIT methods. Moreover, there was a cold front from the north and the wind speed was high on December 24 ( Figure S6 and Figure S7) which makes it more complicated to retrieve BLH from both TBS and lidar data. The example on December 24 suggested that a large difference may be exhibited in BLHs from lidar and TBS data; improving the retrieval method of lidar-BLH does not help to yield accurate BLH from lidar observations under such complicated conditions. Both WCT and FIT methods reasonably captured the diurnal evolution of BLH during the observation period at SACOL (Figure 7). The FIT-BLHs exhibited good consistency with the WCT-BLHs but seemed to be slightly higher than WCT-BLHs, especially during the daytime. This overestimation might be attributed to the influence of elevated aerosols in the residual mixed layers. The standard errors of the WCT-BLHs were generally higher than those of the FIT-BLHs, indicating better stability of the FIT method in determining BLH from lidar observations.

Comparison between BLHs Derived from Lidar Observations and TBS Data
The BLHs derived from lidar observations were compared with those derived from TBS data (Figures 8 and 9). Here, the lidar-BLHs were matched with 58 of the 96 tethered balloon launches. The Pearson correlation coefficients between WCT-BLHs and TBS-BLHs and FIT-BLHs and TBS-BLHs were 0.91 and 0.89, respectively, indicating good consistency between BLHs derived from lidar observations and TBS data. Generally, the Pearson coefficient of correlation and coefficient of determination are used to determine the magnitude of the relation between lidar-BLHs and sounding BLHs. The coefficient of determination, square of the correlation coefficient, is interpreted as the percentage of variance in one variable that is predicted or explained by the other [75]. Some studies have reported coefficient of determination between WCT-BLHs and radiosonde-BLHs of 0.89 for 48 clear sky daytime cases in Houston, USA [49], 0.88 for 23 cases between 13:57 and 21:34 UTC in Maryland, USA [29], and 0.68 for a number of cases in Southern Great Plains site in Oklahoma, USA (36º36' N, 97º29' W) [40]. Moreover, reported correlation coefficients between WCT-BLHs (FIT-BLHs) and radiosonde BLHs include 0.96 (0.94) for 41 cases in early morning and late afternoon transition times (08:00 and 20:00 (Beijing time)) in Lanzhou, China [42], 0.86 (0.88) for 18 cases at 08:00 and 20:00 (Beijing time) in Lanzhou, China [55], and 0.84 for ~65 cases between 06:30 and 18:30 (local time) at Southern Great Plains in Oklahoma, USA [38].

Comparison between BLHs Derived from Lidar Observations and TBS Data
The BLHs derived from lidar observations were compared with those derived from TBS data (Figures 8 and 9). Here, the lidar-BLHs were matched with 58 of the 96 tethered balloon launches. The Pearson correlation coefficients between WCT-BLHs and TBS-BLHs and FIT-BLHs and TBS-BLHs were 0.91 and 0.89, respectively, indicating good consistency between BLHs derived from lidar observations and TBS data. Generally, the Pearson coefficient of correlation and coefficient of determination are used to determine the magnitude of the relation between lidar-BLHs and sounding BLHs. The coefficient of determination, square of the correlation coefficient, is interpreted as the percentage of variance in one variable that is predicted or explained by the other [75]. Some studies have reported coefficient of determination between WCT-BLHs and radiosonde-BLHs of 0.89 for 48 clear sky daytime cases in Houston, USA [49], 0.88 for 23 cases between 13:57 and 21:34 UTC in Maryland, USA [29], and 0.68 for a number of cases in Southern Great Plains site in Oklahoma, USA (36 • 36 N, 97 • 29 W) [40]. Moreover, reported correlation coefficients between WCT-BLHs (FIT-BLHs) and radiosonde BLHs include 0.96 (0.94) for 41 cases in early morning and late afternoon transition times (08:00 and 20:00 (Beijing time)) in Lanzhou, China [42], 0.86 (0.88) for 18 cases at 08:00 and 20:00 (Beijing time) in Lanzhou, China [55], and 0.84 for~65 cases between 06:30 and 18:30 (local time) at Southern Great Plains in Oklahoma, USA [38]. Remote Sens. 2021, 13, x FOR PEER REVIEW 11 of 16  High resolution data (i.e., multiple soundings per day) were able to reveal the diurnal variations in deviations between Lidar-BLHs and TBS-BLHs ( Figure 9). The bias and standard deviation were used to evaluate and quantify the deviations [76,77]. Here, the bias is the difference between the means of lidar-BLHs and TBS-BLHs (Equation S1), the standard deviation (Equation S2) is the root-mean-square value of the departures of the individual pair sample differences from the bias [49], and the t test is used to define the statistical significance of the bias (Equation S3). The overall biases between the TBS-BLHs and the WCT-BLHs (37±79 m) and FIT-BLHs (52±89 m) (not statistically significant; p > 0.05) indicated a very good overall performance by the lidar-based methods. The diurnal variation of the deviations is valuable in validating lidar-BLHs because deviation at a certain time of the day might not be applicable to other times ( Figure S8). The average heights of WCT-BLHs and FIT-BLHs were closest to TBS-BLHs with the bias for 65±106 m (p < 0.05) and 32±148 m (not statistically significant; p > 0.05), respectively, when convection boundary layer heights were at peak, at 14:00. There were small biases between WCT-BLHs and TBS-BLHs for -5±50 m (p < 0.05) and 33±100 m (p < 0.05), respectively, and FIT-BLHs were always more overestimated than TBS-BLHs with the biases of 40±89 m (p < 0.05) and 76±91 m (not statistically significant; p > 0.05) at 11:00 and 17:00, respectively. Moreover, lidar-BLHs were overestimated more obviously with the biases between 41±72 m (p < 0.05) and 66±96 m (not statistically significant; p >  High resolution data (i.e., multiple soundings per day) were able to reveal the diurnal variations in deviations between Lidar-BLHs and TBS-BLHs (Figure 9). The bias and standard deviation were used to evaluate and quantify the deviations [76,77]. Here, the bias is the difference between the means of lidar-BLHs and TBS-BLHs (Equation S1), the standard deviation (Equation S2) is the root-mean-square value of the departures of the individual pair sample differences from the bias [49], and the t test is used to define the statistical significance of the bias (Equation S3). The overall biases between the TBS-BLHs and the WCT-BLHs (37±79 m) and FIT-BLHs (52±89 m) (not statistically significant; p > 0.05) indicated a very good overall performance by the lidar-based methods. The diurnal variation of the deviations is valuable in validating lidar-BLHs because deviation at a certain time of the day might not be applicable to other times ( Figure S8). The average heights of WCT-BLHs and FIT-BLHs were closest to TBS-BLHs with the bias for 65±106 m (p < 0.05) and 32±148 m (not statistically significant; p > 0.05), respectively, when convection boundary layer heights were at peak, at 14:00. There were small biases between WCT-BLHs and TBS-BLHs for -5±50 m (p < 0.05) and 33±100 m (p < 0.05), respectively, and FIT-BLHs were always more overestimated than TBS-BLHs with the biases of 40±89 m (p < 0.05) and 76±91 m (not statistically significant; p > 0.05) at 11:00 and 17:00, respectively. Moreover, lidar-BLHs were overestimated more obviously with the biases between 41±72 m (p < 0.05) and 66±96 m (not statistically significant; p > High resolution data (i.e., multiple soundings per day) were able to reveal the diurnal variations in deviations between Lidar-BLHs and TBS-BLHs (Figure 9). The bias and standard deviation were used to evaluate and quantify the deviations [76,77]. Here, the bias is the difference between the means of lidar-BLHs and TBS-BLHs (Equation (S1)), the standard deviation (Equation (S2)) is the root-mean-square value of the departures of the individual pair sample differences from the bias [49], and the t test is used to define the statistical significance of the bias (Equation (S3)). The overall biases between the TBS-BLHs and the WCT-BLHs (37 ± 79 m) and FIT-BLHs (52 ± 89 m) (not statistically significant; p > 0.05) indicated a very good overall performance by the lidar-based methods. The diurnal variation of the deviations is valuable in validating lidar-BLHs because deviation at a certain time of the day might not be applicable to other times ( Figure S8). The average heights of WCT-BLHs and FIT-BLHs were closest to TBS-BLHs with the bias for 65 ± 106 m (p < 0.05) and 32 ± 148 m (not statistically significant; p > 0.05), respectively, when convection boundary layer heights were at peak, at 14:00. There were small biases between WCT-BLHs and TBS-BLHs for −5 ± 50 m (p < 0.05) and 33 ± 100 m (p < 0.05), respectively, and FIT-BLHs were always more overestimated than TBS-BLHs with the biases of 40 ± 89 m (p < 0.05) and 76 ± 91 m (not statistically significant; p > 0.05) at 11:00 and 17:00, respectively. Moreover, lidar-BLHs were overestimated more obviously with the biases between 41 ± 72 m (p < 0.05) and 66 ± 96 m (not statistically significant; p > 0.05), during the transition period and nighttime (at 02:00, 08:00, and 20:00). Overall, the lidar-based methods slightly overestimated the BLHs most of the time of the day, and the nighttime. In addition, the root-mean-square error (Equation (S4)) between the TBS-BLHs and WCT-BLHs (FIT-BLHs) also showed diurnal variations of 82, 132, 47, 118, 101, and 119 m (82, 97, 92, 142, 116, and 129 m) at 02:00, 08:00, 11:00, 14:00, 17:00, and 20:00, respectively. The overall root-mean-square error of the WCT-BLHs and FIT-BLHs were 105 m and 110 m, respectively. The lidar-BLHs during daytime (i.e., 11:00, 14:00, and 17:00 Beijing time) exhibited very good consistency with the TBS-BLHs, while those during transition times (08:00 and 20:00 Beijing time) and nighttime (02:00 Beijing time) showed more uncertainties when considering both bias and relative errors between lidar-BLHs and TBS-BLHs. Thus, validating lidar-BLHs using observations at transition times may introduce extra errors.
Other studies have reported mean biases between WCT-BLHs and radiosonde-BLHs of 51.1 m for 48 clear sky daytime cases in Houston, USA [49], 280 m for early morning and late afternoon transition times (08:00 and 20:00, Beijing time) in Lanzhou, China [42], and 288 m for 67 daytime cases and 561 m for 64 nighttime cases in Seoul, South Korea [78]. Reported mean biases between FIT-BLHs and radiosonde BLHs include 220 m for early morning and late afternoon transition times (08:00 and 20:00, Beijing time) in Lanzhou, China [42]. The mean overall root-mean-square error of WCT in our study was much lower than that of 105 m in Liu [54] (used eighteen sounding profiles on a single day to validate BLH derived from a ceilometer, the signal of which is like that of a lidar), few studies have used more than six soundings per day to validate the diurnal variation of lidar-BLHs. Overall, our results not only reveal generally consistent results but also showed evident diurnal variations in bias and root-mean-square error for lidar-BLHs and TBS-BLHs. To our best knowledge, our study presents the first comprehensive high-resolution results (3 h time resolution during daytime) on the comparison between lidar-BLHs and TBS-BLHs in China.

Summary and Conclusions
As a vital parameter in atmospheric air pollution studies and predications, highresolution BLHs are usually retrieved from lidar observations. To obtain high-precision and -resolution BLHs and validate their retrieval from lidar observations, simultaneous tethered balloon sounding and lidar observations were conducted at SACOL in Northwest China during 2019-2020 winter with the following main findings.
The Holzworth method and temperature inversion method were selected to derive BLHs (TBS-BLHs) in daytime and nighttime, respectively. The development of the atmospheric boundary was clearly revealed by the high-resolution tethered balloon sounding data. Based on 96 sounding profiles, BLHs were determined as 170, 210, 393, 676, 423, and 190 m at 02:00, 08:00, 11:00, 14:00, 17:00, and 20:00, Beijing time, respectively. These high-resolution TBS-BLHs are valuable for quantifying the dispersion of air pollution in the study region.
The diurnal evolution of BLHs was reasonably captured by both the lidar observationbased WCT and FIT methods, which exhibited correlation coefficients of 0.91 and 0.89, respectively, with TBS-BLHs. Overall, the lidar-based methods slightly overestimated BLHs, which was generally acceptable when considering both the absolute and relative errors compared with the TBS-BLHs. The FIT method exhibited better stability than the WCT method. The lidar-BLHs during daytime exhibited very good consistency with the TBS-BLHs, while those during transition times and nighttime showed more uncertainties when considering both bias and relative errors between lidar-BLHs and TBS-BLHs. Thus, validating lidar-BLHs using observations at transition times may introduce extra errors.
Our results revealed high-resolution and -accuracy diurnal variations of BLHs based on tethered balloon sounding observations at the site of SACOL in Northwest China. The diurnal variation in retrieval errors in lidar-BLHs suggests the importance of validating lidar-based BLHs using simultaneous sounding data. To our best knowledge, our study presents the first comprehensive high-resolution results (3 h time resolution during daytime) on the comparison between lidar-BLHs and TBS-BLHs in China.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/rs13091781/s1, Figure S1: Comparisons between BLHs derived from bulk Richardson number (Rb) method with us=0 and us =0 in daytime, Figure S2: Comparisons between BLHs derived from bulk Richardson number (Rb) method and tethered balloon sounding (TBS) method (a) Holzworth method in daytime and (b) Temperature inversion method in nighttime, Figure S3: Comparisons between BLHs derived from Temperature inversion method and wind shear method in nighttime, Figure S4: Profiles of (a) 02:00 (Beijing time) and (b) 08:00 (Beijing time), 24 December 2019 explain the low stable boundary layer height when the us and vs were not set to zero in the bulk Richardson number (Rb) method in nighttime, Figure Table S1: Detailed information of the tethered sounding system KZXLT-II, Table S2