A Comparison of Atmospheric Boundary Layer Height Determination Methods Using GNSS Radio Occultation Data

: The accurate determination of Atmospheric Boundary Layer Height (ABLH) is crucial in various atmospheric studies and practical applications. In this study, we present a comprehensive comparative analysis of ﬁve distinct methods for estimating ABLH using Global Navigation Satellite System (GNSS) Radio Occultation (RO) data. These methods encompass the use of bending angle and refractivity proﬁles, namely Minimum Gradient methods of the Bending Angle (MGBA) and Refractivity (MGR) proﬁles, breaking point, Wavelet Covariance Transform (WCT), and Double-Parameter Model Function (DPMF). GNSS-RO data from COSMIC-2 and Spire are used. To establish robust validation, radiosonde data are employed as a reference, ensuring the reliability of our ﬁndings. The results reveal notable variations in the performances of these ABLH estimation methods. Speciﬁcally, the MGBA, MGR, breaking point, and DPMF methods exhibit strong correlations with the reference data. Conversely, the WCT method displays weaker correlations, higher biases, and elevated root-mean-square-errors, suggesting limitations in capturing the true ABLH. Furthermore, we remove outlier screening to facilitate a comparison of the differences among the ﬁve methods. The WCT and DPMF methods can detect strong variations in the proﬁles near the Earth’s surface and consider them as ABLH. However, these variations are caused by errors. The MGBA method emerges as a reliable and stable option, while the WCT and DPMF methods should be used with caution due to the lower quality of the GNSS-RO proﬁles near the Earth’s surface.


Introduction
The Atmospheric Boundary Layer (ABL), the lowest layer in the Earth's atmosphere, typically extends upward from the Earth's surface to a few kilometers [1].It plays a crucial role in the vertical exchange of energy, moisture, and pollutants between the Earth's surface and the free troposphere.The Atmospheric Boundary Layer Height (ABLH), which varies in response to temporal and spatial changes, is a critical variable for delineating the characteristics of the ABL within the lower troposphere [2].The ABLH is usually determined by analyzing the vertical profiles of meteorological properties such as temperature, humidity, wind speed, aerosols, or refractivity observed with radiosondes, radars, sodars, and lidars.The Global Navigation Satellite Systems (GNSS) Radio Occultation (RO) technique has significantly advanced the detection and observation of the ABLH by complementing traditional instruments due to its high vertical resolution and all-weather and horizontal global coverage, which advance and complement the capabilities of traditional instruments and enhance our understanding of the ABLH [3][4][5].
Methods for determining ABLH with GNSS-RO observations have been discussed in many previous studies.For example, von Engeln et al. (2005) proposed determining the ABLH from RO using the cut-off altitude defined in the full spectrum inversion method during the processing of Phase-Locked Loop (PLL) data [6].However, Sokolovskiy et al. (2006) [7] pointed out that this method is not reliable because, in many cases, the cut-off altitude determined with PLL tracking data is simply related to the epoch at which the GNSS receiver loses its lock of the signal or starts tracking with large errors due to a low signal-to-noise ratio.Currently, ABLH detection predominantly relies on refractivity and bending angle observations from GNSS-RO data.Several reliable methods have been developed, including the Breaking Point (BP) method, Minimum (Maximum) Gradient (MG) method, Wavelet Covariance Transform (WCT), and regularization method.Sokolovskiy et al. (2006) and Guo et al. (2011) proposed a method for identifying the sharp top of the ABLH using the BP method in the refractivity profile obtained from RO open-loop data and validated the performance of this method by comparing it with high-resolution radiosonde data and ECMWF [7,8].Based on the BP method described by Guo et al. (2011) [8], Chan et al. (2013) [9] implemented stringent constraints on the BP method to effectively handle noise in the refractivity profile and only estimated the sharptopped ABLH, resulting in an ABLH detection frequency of only 31%.The MG method using different atmospheric parameters has been the subject of various research studies.Basha et al. (2009) [10] detected the ABLH from the gradient of the potential temperature, virtual potential temperature, water vapor (mixing ratio), and refractivity profiles.Ao et al. (2012) [11] also compared the ABLH using the minimum gradients of the RO refractivity and specific humidity profiles.They proposed the minimum gradient of the RO-inverted refractivity profile for detecting the ABLH [10,11].Sokolovskiy et al. (2007) provided the MG method using the bending angle profile with less noise than the refractivity profile for determining the ABLH [12].Xie et al. (2012) measured the ABLH over the southeast Pacific Ocean, computing the maximal bending angle lapse instead of the minimum of the refractivity gradient [13].The WCT method, a method for detecting signal mutations, was introduced by Ratnam et al. (2010) to identify the ABLH in refractivity profiles [14].Xu et al. (2023) made advancements to the WCT algorithm by incorporating the available RO-derived ABLH as auxiliary information to detect the ABLH more objectively and focused on investigating the temporal-spatial variation in the Marine Boundary Layer Height (MBLH) in tropics and subtropics [15].The regularization method is a more complex approach that requires higher computational resources and more time compared to other ABLH determination methods, including the L-curve and the Double-Parameter Model Function (DPMF) methods.Yan et al. (2019) utilized the L-curve method, a regularization parameter selection method, to calculate the maximal vertical lapse of the bending angle profile from the Constellation Observing System for Meteorology, Ionosphere, and Climate (COSMIC) bending angle profile [16].Zhou et al. (2020) employed the DPMF method and improved the local convergence problem of the DPMF method, obtaining the bending angle lapse to determine the ABLH [17].
As mentioned above, several methods have been proposed to determine the ABLH, including the MG method, BP method, WCT, and DPMF method using the bending angle or refractivity derived from RO.Each of these methods offers different advantages and disadvantages, depending on the circumstances and the data available for analysis.For example, Ho et al. (2015) declared that the BP method is more susceptible to the perturbations at the top of the inversion layer compared to the MG method, resulting in deriving a higher ABLH value than that from MG method [18].Sokolovskiy et al. (2007) [12] and Ao et al. (2012) [11] highlighted the advantageous properties of the bending angle over refractivity as an observational variable for determining the ABLH, as it has a lower bias and is particularly effective at detecting strong gradients in the atmosphere [11,12].Therefore, it is necessary to conduct a comprehensive investigation on the performances of different methods in determining the ABLH, which will provide an insight for researchers to select an optimal method depending on the circumstances and available data.

Data and Methodology
2.1.COSMIC-2 and Spire COSMIC-2, the second constellation of the program COSMIC, plays an important role in providing accurate bending angle and refraction profiles for atmospheric studies.However, it should be noted that the coverage of COSMIC 2 is currently limited to low latitudes, approximately within ±45 • [19].To solve the problem of insufficient coverage at high latitudes, we integrated RO data from the Spire constellation, which can provide highly accurate bending angle and refractive data with global coverage [20].
The COSMIC-2 and Spire RO data involved in this study are available to freely download from the CDAAC website (https://www.cosmic.ucar.edu/accessed on 24 October 2023).For the purposes of this study, the analysis mainly adopted COSMIC-2 and Spire RO data spanning from 1 March 2022 to 28 February 2023.

IGRA2
The Integrated Global Radiosonde Archive version 2 (IGRA2) is an extensive dataset comprising high-quality sounding observations from over 2800 stations worldwide, freely accessible via the website https://www.ncei.noaa.gov/data/integrated-global-radiosondearchive/access/derived-por/(accessed on 24 October 2023).In this study, the ABLH determined by applying the MG method to the temperature and relative humidity from Radiosonde Observations (RAOB) was used as a reference for validating the RO-inverted ABLH.This validation process was conducted using a spatial window of 300 km and a temporal window of 3 h.

Methods for Deriving ABLH
Both the MG and BP methods are utilized by the COSMIC Data Analysis and Archive Center (CDAAC) to acquire the ABLH from RO profiles.In this study, the ABLHs derived using the MG method with bending angle, MG method with refractivity, BP, WCT, and DPMF methods are denoted as z MGBA , z MGR , z BP , z WCT , and z DPMF .In the aforementioned methods, only the RO data that penetrate below 0.5 km are considered for estimation, and the RO bending angle or refractivity profiles are initially resampled at constant intervals of 5 m using cubic spline interpolation.

MG Method
As per the CDAAC's description, a 300 m window is applied for smoothing the bending angle and refractivity profiles.The gradient is obtained by calculating the finite difference representation, expressed as follows, where X i is the bending angle or refractivity at the ith height and Z i is the altitude.The MG methods for the bending angle and refractivity profiles are abbreviated as MGBA and MGR, respectively.

BP Method
The primary objective of the BP method is to obtain the sharp-topped ABLH.It involves performing a linear regression Az + B on the refractivity profile within a sliding window of 300 m.Then, the lapse rate of coefficient A is calculated, and the ABLH is identified as the height where the maximum lapse rate of A occurs.In addition, for practical applications, the ABLH determined with this method must meet the following criteria [8]: These criteria are essential to ensuring the reliability and accuracy of the ABLH estimated using the BP method.

WCT Method
The WCT algorithm to detect signal mutation was developed by Gamage et al. (1993) [21].
The Haar function h z − b a used by WCT is defined as follows: where z represents the altitude and a and b represent the dilation or spatial extent and center of the Haar function h z−b a , respectively.The dilation, a, is currently unknown and represents the depth of the transition zone [22].The specific value for a will be determined in Section 3.1.For a given refractivity profile f (z), the WCT W f is expressed as: where z t and z b are the top and bottom of the given refractivity profile f (z), respectively.The local maximum of W f corresponds to the ABLH.A peak in W f indicates the presence of a distinct step in the refractivity profile f (z) with a coherent scale of a located at z = b.

DPMF Method
The DPMF method, incorporating numerical regularization, involves a relatively complicated principle and requires significant computational resources.In this context, only a brief overview of the process is presented.For a comprehensive understanding, please consult the specific references [17] cited for further details.Assuming that the bending angle profile is a continuously differentiable function α(z), the gradient of bending angle ϕ(z) = α (z) and z i is the arithmetic altitude sequence with a fixed interval of 5 m.The bending angle gradient satisfies the equation: where, To overcome the ill-posed nature of Equation ( 4), the problem is addressed by transforming the solution into an optimization task, where the objective is to minimize the value of the objective functional.This transformation can be formulated as follows: where X denotes a solution that satisfies the minimum value of the objective functional, the dual regularization parameters, α and β, need to be determined, and the matrix L represents the first derivative operator and has a smooth effect on the gradient x:

Outliers Screening
A sharpness parameter is introduced to quantify disparities and facilitate quality control.This parameter follows a definition similar to the one proposed by Ao et al. (2012) [11].

SP =
X m X RMS (11) where X m and X RMS are the maximum or minimum and Root Mean Square (RMS) of the gradient, lapse of A or W f , respectively.An ABLH with an SP greater than 1.25 should be excluded and regarded as an outlier [9].Furthermore, an RO-derived ABLH is cross-validated using another RO-derived ABLH as auxiliary information.It was identified that there are three markedly distinct ABLHs in this RO profile.As a result, this RO profile will be excluded from the analysis due to these significant dissimilarities, where differences exceeding 0.3 km are deemed as noteworthy.
Lastly, there are practical criteria that need to be considered.The RO profiles must penetrate within a depth below 0.5 km from the surface, and the RO-derived ABLH should not exceed 3.5 km [8,18].

Initial Analysis
CDAAC provides ABLHs calculated using the MG method with bending angle, the MG method with refractivity, and the BP methods [23].Before proceeding with the comparison of the five different methods, we analyze the characteristics of the z MGBA , z MGR , and z BP , calculated by CDAAC.As shown in Figure 1a, some of the ABLHs obtained from the COSMIC-2 data penetrating the below 0.5 km, particularly in the oceanic region, surpass the threshold of 3.5 km.Additionally, Figure 1b illustrates that the z BP do not align with the criteria of |A| > 50 km −1 proposed by Guo et al. (2011) [8].Hence, these values obtained from CDAAC are considered to be raw and have not been filtered according to the specific criteria, making them unreliable and inconsistent with actual atmospheric conditions.Therefore, instead of using the ABLH results provided by CDAAC for the comparative study, we determined the values of z MGBA , z MGR , and z BP using the abovementioned formulae and conducted outlier screening based on the criteria described in Section 2.3.
Based on the experiment conducted by Ratnam et al. (2010) [14], different dilations, ranging from 0.2 km to 0.6 km in intervals of 0.2 km, were selected.To make the dilation more representative of most GNSS-RO profiles, a single GNSS-RO profile was not chosen; instead, all the data on 29 February 2022 were selected to compare the effects of different dilations.Furthermore, it will be compared to the ABLHs calculated using other methods with GNSS-RO profiles.To avoid similarities introduced by the same refractivity profiles, bending angle profiles were used instead.Additionally, considering that the MG method was applied to the determination of AHLH with RO data earlier (than DPMF method) and has received more validation, the MGBA method will be employed.Based on the experiment conducted by Ratnam et al. (2010), different dilations, ranging from 0.2 km to 0.6 km in intervals of 0.2 km, were selected.To make the dilation more representative of most GNSS-RO profiles, a single GNSS-RO profile was not chosen; instead, all the data on 29 February 2022 were selected to compare the effects of different dilations.Furthermore, it will be compared to the ABLHs calculated using other methods with GNSS-RO profiles.To avoid similarities introduced by the same refractivity profiles, bending angle profiles were used instead.Additionally, considering that the MG method was applied to the determination of AHLH with RO data earlier (than DPMF method) and has received more validation, the MGBA method will be employed.
In Figure 2, the proportions of biases > 0.2 km between  and  are 54.79%,14.56%, and 15.04%, corresponding to dilations of 0.2 km, 0.4 km, and 0.6 km, respectively.Therefore, the optimal dilation value should be set to 0.4 km.

Comparison of the RO-and RAOB-Derived ABLHs
To minimize the potential impact of landform and physiognomy on the comparison, the comparison was initiated using data in islands and coastal regions.Therefore, 61 globally distributed RAOB stations (Figure 3) on islands or along coastlines were selected to collocate with the data from GNSS-RO.Additionally, due to GNSS-RO data originating from different RO missions, it was necessary to analyze the COSMIC-2 and Spire data independently to prevent any potential systematic errors.In Figure 2, the proportions of biases > 0.2 km between z MGBA and z WCT are 54.79%,14.56%, and 15.04%, corresponding to dilations of 0.2 km, 0.4 km, and 0.6 km, respectively.Therefore, the optimal dilation value should be set to 0.4 km.Based on the experiment conducted by Ratnam et al. (2010), different dilations, ranging from 0.2 km to 0.6 km in intervals of 0.2 km, were selected.To make the dilation more representative of most GNSS-RO profiles, a single GNSS-RO profile was not chosen; instead, all the data on 29 February 2022 were selected to compare the effects of different dilations.Furthermore, it will be compared to the ABLHs calculated using other methods with GNSS-RO profiles.To avoid similarities introduced by the same refractivity profiles, bending angle profiles were used instead.Additionally, considering that the MG method was applied to the determination of AHLH with RO data earlier (than DPMF method) and has received more validation, the MGBA method will be employed.
In Figure 2, the proportions of biases > 0.2 km between  and  are 54.79%,14.56%, and 15.04%, corresponding to dilations of 0.2 km, 0.4 km, and 0.6 km, respectively.Therefore, the optimal dilation value should be set to 0.4 km.

Comparison of the RO-and RAOB-Derived ABLHs
To minimize the potential impact of landform and physiognomy on the comparison, the comparison was initiated using data in islands and coastal regions.Therefore, 61 globally distributed RAOB stations (Figure 3) on islands or along coastlines were selected to collocate with the data from GNSS-RO.Additionally, due to GNSS-RO data originating from different RO missions, it was necessary to analyze the COSMIC-2 and Spire data independently to prevent any potential systematic errors.

Comparison of the RO-and RAOB-Derived ABLHs
To minimize the potential impact of landform and physiognomy on the comparison, the comparison was initiated using data in islands and coastal regions.Therefore, 61 globally distributed RAOB stations (Figure 3) on islands or along coastlines were selected to collocate with the data from GNSS-RO.Additionally, due to GNSS-RO data originating from different RO missions, it was necessary to analyze the COSMIC-2 and Spire data independently to prevent any potential systematic errors.In this study, using the ABLHs ( ) derived from RAOB atmospheric profiles as a reference, we evaluated the performances of the ABLHs ( ), including  ,  ,  ,  , and  determined from the COSMIC-2 data using five different methods.In this study, using the ABLHs (z RAOB ) derived from RAOB atmospheric profiles as a reference, we evaluated the performances of the ABLHs (z RO ), including z MGBA , z MGR , z BP , z WCT , and z DPMF determined from the COSMIC-2 data using five different methods.As shown in Figure 4, the comparison results indicate a good agreement between the z RAOB and z RO over oceans.The colorbar represents a Gaussian kernel density estimation.Table 1 presents the correlation coefficients between the z RAOB and z RO derived from COSMIC-2, as well as the biases and root-mean-square error (RMSE) of the z RO derived from COSMIC-2 compared to the z RAOB .Higher positive correlation coefficients for z MGBA (0.967), z MGR (0.968), z BP (0.966), z WCT (0.969), and z DPMF (0.967) indicate closer matching between z RO and z RAOB .The differences in the correlation coefficients between z RAOB and z RO were not significant.In terms of bias and RMSE, z RO agreed closely with z RAOB , and were, in descending order: z MGBA , z DPMF , z MGR , z WCT , and z BP .z MGBA exhibited the smallest bias (0.023 km) with an RMSE of 0.173 km and the best agreement with z RAOB , while z BP showed the most noticeable (positive) bias of 0.100 km, with an RMSE of 0.2 km.In this study, using the ABLHs (  ) derived from RAOB atmospheric profiles as reference, we evaluated the performances of the ABLHs (  ), including   ,     ,   , and   determined from the COSMIC-2 data using five different method As shown in Figure 4, the comparison results indicate a good agreement between t   and   over oceans.The colorbar represents a Gaussian kernel density estim tion.Table 1 presents the correlation coefficients between the   and   deriv from COSMIC-2, as well as the biases and root-mean-square error (RMSE) of the   d rived from COSMIC-2 compared to the   .Higher positive correlation coefficients f   (0.967),   (0.968),   (0.966),   (0.969), and   (0.967) indicate clos matching between   and   .The differences in the correlation coefficients betwe   and   were not significant.In terms of bias and RMSE,   agreed closely wi   , and were, in descending order:   ,   ,   ,   , and   .  exhi ited the smallest bias (0.023 km) with an RMSE of 0.173 km and the best agreement wi   , while   showed the most noticeable (positive) bias of 0.100 km, with an RMS of 0.2 km.   Figure 5 shows the scatter density plots of the z RAOB and z RO over oceans determined from the Spire data, suggesting a satisfactory agreement between the z RO and z RAOB , except for the z WCT determined using the WCT method, which had an obvious negative bias compared to the z RAOB .Specifically, some z WCT values were significantly smaller than the values of corresponding z RAOB , with biases of more than ~1.2 km.In Table 2, the correlation coefficient between the z WCT calculated from Spire and z RAOB is 0.686.
Compared to the correlation between the z WCT derived from COSMIC-2 and z RAOB (0.969), the correlation weakened significantly.Moreover, the bias and RMSE were relatively large.In terms of bias and RMSE, z RO agreed closely with z RAOB , and were, in descending order: z MGBA , z DPMF , z MGR , z BP , and z WCT .The bias and RMSE for z MGBA remained the smallest compared to those from COMSIC-2, showing the closest agreement with z RAOB .
mained the smallest compared to those from COMSIC-2, showing the closest agreem with  .Compared to RAOB stations on islands and near coastal regions, there are more RAOB stations on land, and the ABLH over land shows a greater variability than that over the ocean.However, inland RAOB stations are more affected by landform and physiognomy.Therefore, we only considered stations with altitudes below 0.5 km, located deep inland at least 300 km away from the coast, to avoid the influence of coastal and oceanic areas.The distribution of inland RAOB stations is shown in Figure 6.
Figure 7 displays the scatter density plots of the z RAOB and z RO calculated from the COSMIC-2 data for inland regions.Table 3 shows the statistical analysis results of the z RAOB and z RO inland from the COSMIC-2 data.Combining the statistical results from Table 3, the main findings and the results from the z RO over oceans for Spire were similar.The number of RO profiles over oceans was relatively limited, which may explain why smaller z WCT values were not found in the COSMIC-2 oceanic data.This suggests that the cause of the significantly lower z WCT compared to z RAOB was likely inherent to the WCT method itself, rather than the Spire RO data.
Compared to RAOB stations on islands and near coastal regions, there are more RAOB stations on land, and the ABLH over land shows a greater variability than that over the ocean.However, inland RAOB stations are more affected by landform and physiognomy.Therefore, we only considered stations with altitudes below 0.5 km, located deep inland at least 300 km away from the coast, to avoid the influence of coastal and oceanic areas.The distribution of inland RAOB stations is shown in Figure 6.   3 shows the statistical analysis results of the  and  inland from the COSMIC-2 data.Combining the statistical results from Table 3, the main findings and the results from the  over oceans for Spire were similar.The number of RO profiles over oceans was relatively limited, which may explain why smaller  values were not found in the COSMIC-2 oceanic data.This suggests that the cause of the significantly lower  compared to  was likely inherent to the WCT method itself, rather than the Spire RO data.3 shows the statistical analysis results of the  and  inland from the COSMIC-2 data.Combining the statistical results from Table 3, the main findings and the results from the  over oceans for Spire were similar The number of RO profiles over oceans was relatively limited, which may explain why smaller  values were not found in the COSMIC-2 oceanic data.This suggests that the cause of the significantly lower  compared to  was likely inherent to the WCT method itself, rather than the Spire RO data.In Figure 8, the scatter density plots of the z RAOB and z RO inland calculated from the Spire data are shown.Table 4 presents the statistical analysis results of the z RAOB and z RO from the Spire data for inland regions.In Figure 8, the ABLHs calculated by Spire in the range from 1 to 1.5 km are denser compared to COSMIC-2.There is a greater quantity of data in the range from 1 to 1.5 km, in contrast to the inland ABLHs obtained from COSMIC-2.This was mainly because COSMIC-2 is concentrated in mid-latitude regions, whereas Spire offers global coverage.The ABLH in mid-latitude regions was higher compared to that in high-latitude regions [24]. 0.970 0.023 0.157 In Figure 8, the scatter density plots of the  and  inland calculated from the Spire data are shown.Table 4 presents the statistical analysis results of the  and  from the Spire data for inland regions.In Figure 8, the ABLHs calculated by Spire in the range from 1 to 1.5 km are denser compared to COSMIC-2.There is a greater quantity o data in the range from 1 to 1.5 km, in contrast to the inland ABLHs obtained from COS MIC-2.This was mainly because COSMIC-2 is concentrated in mid-latitude regions whereas Spire offers global coverage.The ABLH in mid-latitude regions was higher com pared to that in high-latitude regions [24].Based on the statistical data presented in Tables 1-4, it can be observed that there were no significant differences in the biases and RMSE of the  from Spire and  Based on the statistical data presented in Tables 1-4, it can be observed that there were no significant differences in the biases and RMSE of the z RO from Spire and COSMIC-2 compared to the z RAOB , except in the case of the WCT method.As multiple reports and studies have already compared and evaluated data from both COSMIC-2 and Spire, and considering the high similarity in the bending angle and refractivity accuracy [25,26], we will not discuss these aspects separately in the following content.
Figure 9 illustrates the biases and RMSE of the z RO over oceans in the Northern Hemisphere using the five different methods across four seasons (MAM: from March to May, JJA: from June to August, SON: from September to November, and DJF: from December to February).Specially, z * WCT represents the z WCT without a bias of more than 1.2 km with z RAOB .In all four seasons, z WCT consistently had the largest bias and RMSE among the five methods.Notably, the biases of z WCT and z * WCT exhibited seasonal variability, with the highest biases being observed during JJA and the lowest during DJF.Conversely, the RMSE values for the z MGBA , z MGR , z BP , and z DPMF did not exhibit pronounced seasonal variations.The biases of z * WCT were smaller than those of z WCT , which was due to the removal of data with biases exceeding 1.2 km.The RMSE of z WCT was then relatively close to that of other z RO .However, the biases and RMSE of z * WCT still showed some seasonal characteristics, possibly because larger biases were not completely removed, and the lower number of collocated RO profiles over oceans made it difficult to effectively reduce errors through averaging.It is possible that this was caused by the lower quality of the GNSS-RO refractivity profiles near the Earth's surface.Near the surface, super-refraction and signal tracking errors are more likely to occur, leading to negative biases in refractivity [27][28][29].Environments with high humidity are more prone to inducing super-refraction [30].As a result, the bias and RMSE of z WCT are more likely to exhibit seasonal variations.Further analysis will be conducted in Section 3.4.
isphere using the five different methods across four seasons (MAM: from to May, JJA: from June to August, SON: from September to November, and DJF: from December to February).Specially,  * represents the  without a bias of more than 1.2 km with . In all four seasons,  consistently had the largest bias and RMSE among the five methods.Notably, the biases of  and  * exhibited seasonal variability, with the highest biases being observed during JJA and the lowest during DJF.Conversely, the RMSE values for the  ,  ,  , and  did not exhibit pronounced seasonal variations.The biases of  * were smaller than those of  , which was due to the removal of data with biases exceeding 1.2 km.The RMSE of  was then relatively close to that of other  .However, the biases and RMSE of  * still showed some seasonal characteristics, possibly because larger biases were not completely removed, and the lower number of collocated RO profiles over oceans made it difficult to effectively reduce errors through averaging.It is possible that this was caused by the lower quality of the GNSS-RO refractivity profiles near the Earth's surface.Near the surface, super-refraction and signal tracking errors are more likely to occur, leading to negative biases in refractivity [27][28][29].Environments with high humidity are more prone to inducing super-refraction [30].As a result, the bias and RMSE of  are more likely to exhibit seasonal variations.Further analysis will be conducted in Section 3.4.In Figure 10, the deviation and RMSE between  and  within the Northern Hemisphere's inland regions are presented.It can be observed that the  * in inland areas did not exhibit significant seasonal variations, due to a larger number of profiles, which helped to reduce errors through averaging.From Tables 1-4 and Figures 9 and 10, it is evident that  performed the best and was closest to the RAOB measurements.In Figure 10, the deviation and RMSE between z RO and z RAOB within the Northern Hemisphere's inland regions are presented.It can be observed that the z * WCT in inland areas did not exhibit significant seasonal variations, due to a larger number of profiles, which helped to reduce errors through averaging.From Tables 1-4 and Figures 9 and 10, it is evident that z MGBA performed the best and was closest to the RAOB measurements.

Comparison of the RO-Derived ABLHs
The comparisons between the ABLHs derived from the RO data and the RAOB Section 3.1 are restricted only to the vicinity of RAOB stations.To further compare th ALBHs derived from RO data on a large spatial scale, we conducted a comprehensi investigation of MBLHs globally.Based on the above analysis,  exhibited the sma est error.Therefore, we used  as a reference value to analyze other  .Figure 11 shows the annual-averaged differences of  ,  ,  , an  w.r.t. over oceans for the period from 1 March 2022 to 28 February 2023 in given latitude-longitude bin with a resolution of 1° × 1°.In ascending order w.r.t. the d ferences from  , they were  ,  ,  , and  .Specifically, the differenc

Comparison of the RO-Derived ABLHs
The comparisons between the ABLHs derived from the RO data and the RAOB in Section 3.1 are restricted only to the vicinity of RAOB stations.To further compare the ALBHs derived from RO data on a large spatial scale, we conducted a comprehensive investigation of MBLHs globally.Based on the above analysis, z MGBA exhibited the smallest error.Therefore, we used z MGBA as a reference value to analyze other z RO .
Figure 11 shows the annual-averaged differences of z MGR , z BP , z WCT , and z DPMF w.r.t.z MGBA over oceans for the period from 1 March 2022 to 28 February 2023 in a given latitudelongitude bin with a resolution of 1 • × 1 • .In ascending order w.r.t. the differences from z MGBA , they were z WCT , z MGR , z DPMF , and z BP .Specifically, the differences z WCT and z MGR were negative, while those for z DPMF and z BP were positive.The z WCT values were relatively lower than the other z RO values.This was because z WCT exhibited a significant negative bias, as presented in Section 3.2.
The data are grouped into four different seasons (MAM: from March to May, JJA: from June to August, SON: from September to November, and DJF: from December to February). * represents the  without a bias of more than 1.2 km with  .

Comparison of the RO-Derived ABLHs
The comparisons between the ABLHs derived from the RO data and the RAOB in Section 3.1 are restricted only to the vicinity of RAOB stations.To further compare the ALBHs derived from RO data on a large spatial scale, we conducted a comprehensive investigation of MBLHs globally.Based on the above analysis,  exhibited the smallest error.Therefore, we used  as a reference value to analyze other  .Figure 11 shows the annual-averaged differences of  ,  ,  , and  w.r.t. over oceans for the period from 1 March 2022 to 28 February 2023 in a given latitude-longitude bin with a resolution of 1° × 1°.In ascending order w.r.t. the differences from  , they were  ,  ,  , and  .Specifically, the differences for  and  were negative, while those for  and  were positive.The  values were relatively lower than the other  values.This was because  exhibited a significant negative bias, as presented in Section 3.2.The boxplot in Figure 12 graphically shows the distribution of the MBLHs determined by the five methods from 1 March 2022 to 28 February 2023.The box extends from the first quartile (Q1) to the third quartile (Q3), and the yellow line indicates the median value.More detailed information on the median, Q1, and Q3 values for the MBLHs obtained using each method can be found in Table 5.Among the five methods, the distributions of z MGBA and z DPMF showed a better agreement with each other than with the others, with a median value of about 1.6 km.Compared to z MGBA and z DPMF , the median value of z BP (1.820) was relatively high, while the median value of z MGR (1.598) wazs relatively low.The values of z WCT , on the other hand, were significantly lower than those of the other methods, which is consistent with the results presented in Section 3.2.

Analysis on the Biases of the RO-derived ABLHs
The experimental analysis mentioned above involved the outlier screening of  using the five methods. was consistently lower than the ABLHs obtained using other methods.This discrepancy may be attributed to the sensitivity of the WCT method to refractivity mutations near the Earth's surface.In this case, outlier screening was re moved to facilitate the comparison of the differences among the five methods, especially the WCT method.
Compared to Figure 12 and Table 5, as shown in Figure 13 and Table 6, the median values of the MBLH calculated using the five methods decreased.This was due to th outlier screening process, which removed data with larger differences.MBLH values nea the Earth's surface are more numerous and contain more errors, making them more likely to be removed.Without outlier screening, the number of smaller values for  would

Analysis on the Biases of the RO-Derived ABLHs
The experimental analysis mentioned above involved the outlier screening of z RO using the five methods.z WCT was consistently lower than the ABLHs obtained using other methods.This discrepancy may be attributed to the sensitivity of the WCT method to refractivity mutations near the Earth's surface.In this case, outlier screening was removed to facilitate the comparison of the differences among the five methods, especially the WCT method.
Compared to Figure 12 and Table 5, as shown in Figure 13 and Table 6, the median values of the MBLH calculated using the five methods decreased.This was due to the outlier screening process, which removed data with larger differences.MBLH values near the Earth's surface are more numerous and contain more errors, making them more likely to be removed.Without outlier screening, the number of smaller values for z WCT would increase noticeably, followed by z DPMF .In contrast, the boxplots of the other three methods show relatively small variations after the outlier screening (as seen in Figure 12).The DPMF and WCT methods attempt to utilize data near the lower troposphere more extensively than the other three methods.However, the RO profiles near the lower troposphere are associated with larger errors and significant variations, resulting in the unreliability of the z WCT and z DPMF .
Atmosphere 2023, 14, x FOR PEER REVIEW 14 of 18 increase noticeably, followed by  .In contrast, the boxplots of the other three methods show relatively small variations after the outlier screening (as seen in Figure 12).The DPMF and WCT methods attempt to utilize data near the lower troposphere more extensively than the other three methods.However, the RO profiles near the lower troposphere are associated with larger errors and significant variations, resulting in the unreliability of the  and  .As shown in Figure 14, a GNSS-RO event profile of bending angle (red line) and refractivity (blue line) was selected to further analyze the biases of the RO-derived ABLHs.The WCT and DPMF methods detected strong variations in the profiles near the Earth's surface and considered them as ABLHs.Through comparison to z RAOB , it became evident that these pronounced fluctuations were more likely caused by errors.Using other three z RO as auxiliary information, the heights of the minimum points closest to the values of other three z RO were considered as the ABLHs obtained using WCT and DPMF (corresponding to z WCT and z DPMF ).
Atmosphere 2023, 14, x FOR PEER REVIEW 1 Basing their values on the other three  , the occurrence of significantly low ues in the case of the  and  was considered as an outlier.We conducted tistical analysis on the occurrence of lower values for  and  from 1 Marc to 28 February 2023, and found that the probabilities were 10.50% and 4.24%, respec Basing their values on the other three z RO , the occurrence of significantly lower values in the case of the z WCT and z DPMF was considered as an outlier.We conducted a statistical analysis on the occurrence of lower values for z WCT and z DPMF from 1 March 2022 to 28 February 2023, and found that the probabilities were 10.50% and 4.24%, respectively.

Discussion
Overall, our study focused on the comparison and evaluation of ABLH estimation methods using data from various sources, including GNSS-RO and RAOB.We analyzed five different methods: MGBA, MGR, BP, WCT, and DPMF.
First, we analyzed the raw ABLHs provided by the CDAAC and identified the significant discrepancies and inconsistencies in the z MGBA , z MGR , and z BP values that did not align with the established criteria.As a result, we opted to derive z MGBA , z MGR , and z BP using our own calculations and applied outlier screening based on predefined criteria.This ensured the reliability and accuracy of the ABLH data used in our comparative study.Furthermore, an analysis of the optimal dilation value for WCT was conducted.Based on extensive data statistics, the optimal value was determined to be 0.4 km.
In terms of correlation, z MGBA , z MGR , z BP , and z DPMF demonstrated strong positive correlations with z RAOB .When considering biases and RMSE, z MGBA exhibited minimal errors, showing relatively small biases and RMSE values that closely aligned with the reference z RAOB .In contrast, the WCT method stood out with a little weaker correlation, noticeable bias, and higher RMSE.z WCT occasionally exhibited some low values, causing a decrease in correlation, as well as increases in bias and RMSE.
A seasonal analysis revealed that the bias in the ABLH estimations for z MGBA , z MGR , z BP , and z DPMF exhibited no obvious seasonal variability.Notably, the biases of z WCT exhibited seasonal variability, with the highest biases being observed during JJA and the lowest during DJF.It is possible that this was caused by the lower quality of the GNSS-RO refractivity profiles near the Earth's surface.Near the surface, super-refraction and signal tracking errors are more likely to occur, leading to negative biases in refractivity [27][28][29].Environments with high humidity are more prone to inducing super-refraction [30].As a result, the bias and RMSE of z WCT were more likely to exhibit seasonal variations.
Furthermore, we investigated the global distribution of the MBLHs obtained using the five methods.While z MGBA ,z DPMF and z MGR exhibited minimal differences, the z WCT values were relatively lower than the other z RO values, while the z BP values were relatively higher than the other z RO values.
Eventually, we removed the outlier screening to facilitate a comparison of the differences among the five methods.The WCT and DPMF methods could detect strong variations in the profiles near the Earth's surface and considered them as ABLH.Through a comparison to z RAOB , it became evident that these pronounced fluctuations were more likely caused by errors.
In summary, these findings have important implications for selecting appropriate methods for ABLH estimation.The MGBA method emerged as a reliable and stable option, while the WCT and DPMF methods should be used with caution due to the lower quality of the GNSS-RO refractivity profiles near the Earth's surface.Furthermore, the DPMF method has its shortcomings, including lengthy computations and the possibility of failed iterative calculations.The MGR and BP methods also showed promising performances.

Conclusions
In this work, our study compared five ABLH estimation methods using GNSS-RO and RAOB data.Based on the above assessment and analysis, the following conclusions can be drawn:

•
z MGBA exhibited minimal errors, showing relatively small biases and RMSE values that closely aligned with the reference z RAOB .

•
z WCT occasionally exhibited some low values, causing a decrease in correlation, as well as increases in bias and RMSE.

•
A seasonal analysis revealed that the bias in the ABLH estimations for z MGBA , z MGR , z BP , and z DPMF exhibited no obvious seasonal variability.Notably, the biases of z WCT exhibited seasonal variability, with the highest biases being observed during JJA and the lowest during DJF.

•
The WCT and DPMF methods could detect strong variations in the profiles near the Earth's surface and considered them as ABLH.However, these variations were caused by errors.They should be used with caution due to the lower quality of the GNSS-RO bending angle or refractivity profiles near the Earth's surface.
This research comprehensively compared ABLH estimation methods, with MGBA standing out for its exceptional performance.It offers highly accurate ABLH estimates that closely agree with RAOB-derived ABLHs, providing robust support for precise atmospheric boundary layer height estimation.

Figure 3 .
Figure 3. Location mapping: selected stations on islands or along coastlines from IGRA2 dataset.

Figure 3 .
Figure 3. Location mapping: selected stations on islands or along coastlines from IGRA2 dataset.

Figure 3 .
Figure 3. Location mapping: selected stations on islands or along coastlines from IGRA2 dataset.

Figure 4 .
Figure 4. Scatter density plots between the z RAOB and z RO over oceans for COSMIC-2 with colorbar representing the Gaussian kernel density estimation: (a) z MGBA , (b) z MGR , (c) z BP , (d) z WCT , and (e) z DPMF .

Figure 5 .
Figure 5. Scatter density plots between the  and  over oceans for Spire with colorbar resenting the Gaussian kernel density estimation: (a)  , (b)  , (c)  , (d)  , and  .

Figure 5 .
Figure 5. Scatter density plots between the z RAOB and z RO over oceans for Spire with colorbar representing the Gaussian kernel density estimation: (a) z MGBA , (b) z MGR , (c) z BP , (d) z WCT , and (e) z DPMF .

Figure 6 .
Figure 6.Location mapping: selected stations over land from IGRA2 dataset.

Figure 7
Figure 7 displays the scatter density plots of the  and  calculated from the COSMIC-2 data for inland regions.Table3shows the statistical analysis results of the  and  inland from the COSMIC-2 data.Combining the statistical results from Table3, the main findings and the results from the  over oceans for Spire were similar.The number of RO profiles over oceans was relatively limited, which may explain why smaller  values were not found in the COSMIC-2 oceanic data.This suggests that the cause of the significantly lower  compared to  was likely inherent to the WCT method itself, rather than the Spire RO data.

Figure 6 .
Figure 6.Location mapping: selected stations over land from IGRA2 dataset.

Figure 6 .
Figure 6.Location mapping: selected stations over land from IGRA2 dataset.

Figure 7
Figure 7 displays the scatter density plots of the  and  calculated from the COSMIC-2 data for inland regions.Table3shows the statistical analysis results of the  and  inland from the COSMIC-2 data.Combining the statistical results from Table3, the main findings and the results from the  over oceans for Spire were similar The number of RO profiles over oceans was relatively limited, which may explain why smaller  values were not found in the COSMIC-2 oceanic data.This suggests that the cause of the significantly lower  compared to  was likely inherent to the WCT method itself, rather than the Spire RO data.

Figure 7 .
Figure 7. Scatter density plots between the z RAOB and z RO inland for COSMIC-2 with colorbar representing the Gaussian kernel density estimation: (a) z MGBA , (b) z MGR , (c) z BP , (d) z WCT , and (e) z DPMF .

Figure 8 .
Figure 8. Scatter density plots between the z RAOB and z RO inland for Spire with colorbar representing the Gaussian kernel density estimation: (a) z MGBA , (b) z MGR , (c) z BP , (d) z WCT , and (e) z DPMF .

Figure 9 .
Figure 9. Biases (a) and RMSE (b) of  over oceans in the Northern Hemisphere using five methods.The data are grouped into four different seasons (MAM: from March to May, JJA: from June to August, SON: from September to November, and DJF: from December to February). * represents the  without a bias of more than 1.2 km with  .

Figure 9 .
Figure 9. Biases (a) and RMSE (b) of z RO over oceans in the Northern Hemisphere using five methods.The data are grouped into four different seasons (MAM: from March to May, JJA: from June to August, SON: from September to November, and DJF: from December to February).z * WCT represents the z WCT without a bias of more than 1.2 km with z RAOB .

Atmosphere 2023 ,
14, x FOR PEER REVIEW 12 of

Figure 10 .
Figure 10.Biases (a) and RMSE (b) of  inland in the Northern Hemisphere using five method The data are grouped into four different seasons (MAM: from March to May, JJA: from June August, SON: from September to November, and DJF: from December to February). * repr sents the  without a bias of more than 1.2 km with  .

Figure 10 .
Figure 10.Biases (a) and RMSE (b) of z RO inland in the Northern Hemisphere using five methods.The data are grouped into four different seasons (MAM: from March to May, JJA: from June to August, SON: from September to November, and DJF: from December to February).z * WCT represents the z WCT without a bias of more than 1.2 km with z RAOB .

Figure 12 .
Figure 12.The boxplot of the MBLHs using the five methods from 1 March 2022 to 28 February 2023

Figure 12 .
Figure 12.The boxplot of the MBLHs using the five methods from 1 March 2022 to 28 February 2023.

Figure 13 .Figure 13 .
Figure 13.The boxplot of the MBLHs without outlier screening using the five methods from 1 March 2022 to 28 February 2023.As shown in Figure 14, a GNSS-RO event profile of bending angle (red line) and refractivity (blue line) was selected to further analyze the biases of the RO-derived ABLHs.The WCT and DPMF methods detected strong variations in the profiles near the Earth's surface and considered them as ABLHs.Through comparison to  , it became evident that these pronounced fluctuations were more likely caused by errors.Using other three  as auxiliary information, the heights of the minimum points closest to the values of Figure 13.The boxplot of the MBLHs without outlier screening using the five methods from 1 March 2022 to 28 February 2023.

Figure 14 .
Figure 14.The GNSS-RO event profile of bending angle (red line) and refractivity (blue lin five  . and  represented using  with the other three methods as auxiliary mation to determine the correct ABLH by selecting the minimum point heights from the three ods closest to them.

Figure 14 .
Figure 14.The GNSS-RO event profile of bending angle (red line) and refractivity (blue line) with five z RO .z WCT and z DPMF represented using z RO with the other three methods as auxiliary information to determine the correct ABLH by selecting the minimum point heights from the three methods closest to them.

Table 1 .
The statistical analysis of the data z RAOB and z RO over oceans derived from COSMIC-2.

Table 2 .
The statistical analysis of the data  and  over oceans derived from Spire.

Table 2 .
The statistical analysis of the data z RAOB and z RO over oceans derived from Spire.

Table 3 .
The statistical analysis of the data z RAOB and z RO inland derived from COSMIC-2.

Table 4 .
The statistical analysis of the data  and  inland derived from Spire

Table 4 .
The statistical analysis of the data z RAOB and z RO inland derived from Spire.