Correcting negatively biased refractivity below ducts in GNSS radio occultation: an optimal estimation approach towards improving planetary boundary layer (PBL) characterization

Global Navigation Satellite System (GNSS) radio occultation (RO) measurements are promising in sensing the vertical structure of the Earth’s planetary boundary layer (PBL). However, large refractivity changes near the top of PBL can cause ducting and lead to a negative bias in the retrieved refractivity within the PBL (below ∼ 2 km). To remove the bias, a reconstruction method with assumption of linear structure inside the ducting layer models has been proposed by Xie et al. (2006). While the negative bias can be reduced drastically as demonstrated in the simulation, the lack of high-quality surface refractivity constraint makes its application to real RO data difficult. In this paper, we use the widely available precipitable water (PW) satellite observation as the external constraint for the bias correction. A new framework is proposed to incorporate optimization into the RO reconstruction retrievals in the presence of ducting conditions. The new method uses optimal estimation to select the best refractivity solution whose PW and PBL height best match the externally retrieved PW and the known a priori states, respectively. The near-coincident PW retrievals from AMSR-E microwave radiometer instruments are used as an external observational constraint. This new reconstruction method is tested on both the simulated GNSS-RO profiles and the actual GNSS-RO data. Our results show that the proposed method can greatly reduce the negative refractivity bias when compared to the traditional Abel inversion.


Introduction
The planetary boundary layer (PBL) is the lowest layer of the atmosphere (∼ 2 km) and couples the surface to the free troposphere.Influenced mainly by surface friction, solar radiation, and turbulent transport of moisture, the PBL controls the energy distribution from the surface into the atmosphere.Through the turbulent winds along with cumuliform and stratiform clouds formation, the PBL can greatly affect the local weather as well as the global climate (Garratt, 1992).Due to its importance to the weather prediction community, the PBL has been extensively studied with various sounding techniques for several decades.
Among the existing probing technologies, Global Navigation Satellite System (GNSS) radio occultation (RO), which provides high-resolution atmospheric vertical profiles, has been used to characterize the PBL in recent years (von Engeln et al., 2005;Sokolovskiy et al., 2007).GNSS RO is a limb-sounding technique that precisely measures the GNSS signal phase delay received by low-Earth-orbiting satellites, through which the bending angle and accurate atmospheric refractivity profiles can be retrieved (Kursinski et al., 1997).As a remote sensing technique, GNSS RO acquires measurements over remote regions including marine areas where the PBL can be crucial for weather (Beljaars and Viterbo, 1998) and climate modeling (Zeng et al., 2004) and cannot be achieved by traditional radiosonde, tower, and field observations.A high vertical resolution of approximately 100 m (Gorbunov et al., 2004) is another advantage of GNSS Published by Copernicus Publications on behalf of the European Geosciences Union.
RO over other passive remote sensing instruments (Curran, 1989).Additionally, the L-band GNSS RO signals can penetrate through clouds and precipitation (Solheim et al., 1999), which are common at the height at the top of the PBL.These features make GNSS RO a valuable tool for sensing the PBL (Guo et al., 2011;Ao et al., 2012;Xie et al., 2012;Chan and Wood, 2013;Ho et al., 2015).
However, it is known that the large refractivity change associated with a strong inversion layer at the top of the PBL can cause severe negative biases in RO refractivity measurements (N -bias) (Sokolovskiy, 2003;Xie et al., 2006;Ao, 2007).The large temperature and moisture changes near the top of the PBL could lead to a sharp negative refractivity gradient such that the radius of curvature of the signal path can become less than the radius of the Earth.This phenomenon, called ducting, occurs when the refractivity gradient dN/dr −157 (N -units km −1 ) and can be frequently observed in the subtropics below 2 km.The layer where the ducting occurred is called the ducting layer.Due to the transmitter and receiver geometry of GNSS RO, the tangent point of each ray path is never located within the ducting layers, which theoretically can "trap" the signal whose tangent point is inside.As a result, the GNSS-RO bending angle measurements will lose the information inside the ducting layer, in which the information cannot be recovered using solely GNSS-RO observations.The standard Abel inversion of the bending angle profile will always lead to a profile with no ducts and can cause a negative N-bias as large as 15 % below the ducting layer (Xie et al., 2010).Correcting the Nbias within the PBL is essential towards the use of RO in studying the vertical structure within the PBL.While the weather analyses can assimilate RO bending angles, which are unaffected by the refractivity bias caused by ducting, it is not clear that the analyses can optimally handle these highvertical-resolution measurements.In addition, the analyses may be strongly affected by bias in the model, as evidenced by the low PBL height over the stratocumulus regions (Xie et al., 2012).Therefore, it is of great scientific interests to retrieve an unbiased PBL refractivity based on observations only.
To mitigate the N-bias and reconstruct refractivity profiles inside the boundary layer, a reconstruction method was proposed by Xie et al. (2006), hereinafter referred to as (Xie06).The paper confirmed that an infinite number of refractivity profiles correspond to one bending angle profile in the presence of ducting conditions.A nonlinear function used to describe the continuum of refractivity solutions was derived based on the Abel-retrieved refractivity profile.Choosing the correct parameter and profile from the continuum, however, depends on two assumptions that cannot be easily fulfilled.First, to use the surface refractivity constraint, the RO bending angle measurements are implicitly assumed to cover all altitudes and stop exactly at the Earth's surface.However, the real RO bending angle profiles often do not reach the surface due to a combination of receiver measurement errors and at-mospheric variabilities (Ao et al., 2012).Second, the reconstruction method assumes that the top height of the ducting layer can be determined accurately.However, due to the high variability in the bending angle, identifying the impact parameter of the ducting layer in the real occultation could be challenging.
In this paper we present a new and improved reconstruction method that implements optimal estimation along with external measurements of precipitable water (PW) based on a modification of the Xie06 approach.In Sect.2, the ducting effects and the reconstruction method in Xie06 are reviewed and our new approach using optimal estimation is described.The results for radiosonde observation (RAOB) using the optimal estimation approach and the comparison with different PW observation sources are presented in Sect.3. In Sect. 4 we validate actual GNSS-RO data results using the proposed reconstruction method.The summaries and conclusions are provided in Sect. 5.

N -bias
Under the assumption of a spherically symmetric atmosphere, the impact parameter a of a single ray path can be defined as where n is the refractive index, r is the distance from the center of curvature to each point of the ray path, and φ is the angle between the ray path and the radial vector.The accumulated bending angle of a GNSS-RO ray path can be calculated for a refractivity profile as (Fjeldbo et al., 1971) where α is the bending angle and r t is the radius of the ray path at the tangent point.To retrieve the refractivity information from the bending angle measurement, Eq. ( 2) can be simplified by using the impact parameter defined as a = n (r t ) r t (φ = 90 • in Eq. ( 1) at the tangent point) and assuming that the function x = n (r) r is monotonically increasing with r: so that the refractive index can be derived analytically as (Fjeldbo et al., 1971) which is called the Abel-inversion integral.The Abel inversion is extensively used in RO retrievals, based on the fact that in most cases the one-to-one relationship between the derived refractivity profile and the measured bending angle profile is valid.However, this relationship breaks down when ducting occurs (Sokolovskiy, 2003).

Ducting effects
The refractivity in the neutral atmosphere is related to atmospheric temperature, pressure, and the water vapor pressure with the following equation (Smith and Weintraub, 1953): where N = (n−1)×10 6 is the refractivity in N -units, n is the refractive index, p is pressure in millibar, T is temperature in degrees Kelvin, and e is the water vapor pressure in millibar.Due to the large change in temperature and moisture, the refractivity decreases rapidly across the top of the PBL as seen in Fig. 1a between the height h m and h t .The ducting condition occurs when the refractivity gradient exceeds the critical refraction; i.e., dN/dr − 157 (N-units km −1 ), where x (r) = n (r) r is no longer a monotonic function with respect to r.The height h is defined (Sokolovskiy, 2003) as where r e is the radius of curvature at the Earth's surface.As illustrated in Fig. 1b, the function h(x), shown in black, can be divided into four intervals: h 1 (x), where x increases from the surface to the height h b to reach x b ; h 2 (x), where x further increases to x m , from h b to h m ; h 3 (x), which is the ducting layer with refractivity gradient exceeding critical refraction, where x decreases from x m at h m back to x b at h t ; and h 4 (x), where x increases again from h t and the monotonic relation between x and r is restored.Within the intervals of h 2 (x) and h 3 (x), the changing signs of the slopes result in a non-monotonic relationship between h and x.
For simplicity, here, we define the trapping layer, which includes the ducting layer and the layer underneath from h b to h t where the monotonic characteristic of h(x) vanishes.Inside the trapping layer, the refractivity gradient is large enough to trap the signal with a tangent point between h b and h t and cause an infinite bending angle in its ray path.Because of the geometry of GNSS RO, in which both the transmitter and the receiver are located outside the Earth's atmosphere, the tangent points of the received signals do not appear inside the trapping layer.In other words, the information between h b and h t is lost in the received signal, and the bending angle observation is not able to cover this gap.This gap can be noticed by examining Eq. (2).When evaluating the bending angle with Eq. (2) inside the trapping layer, h b < r t −r e < h t , the term n(r)r above the height r t becomes less than n(r t )r t because of the negative gradient of x between h m and h t .This would lead to a negative value inside the square root in Eq. ( 2) and the solution would be a complex number for the bending angle inside the trapping layer which is unphysical.However, all the bending angles below h b can still be evaluated as the monotonic relationship of h(x) is in place.Namely, the rays with tangent points below the trapping layer can still penetrate through the trapping layer and arrive at the receiver, and so the bending angle of these rays can still be calculated by GNSS RO.
Although the missing bending angle measurements in the trapping layer cause a gap in the α − r relationship, the α − a relationship, will remain seamless because the x value at the top and the bottom of the trapping layer are identical (e.g., x b in Fig. 1b).Therefore, we can still apply the standard Abel inversion upon the bending angle profile in the presence of the ducting.However, since the x(r) function is not monotonic within the trapping layer, the standard Abel inversion Eq. ( 4), which assumes monotonic x, will lead to erroneous refractivity results below h t .While the refractivity retrieval remains valid above h t , the bending contribution inside the trapping layer is missing from the standard Abel inversion below x b .Consequently, the Abel-retrieved refractivity below x b will be negatively biased (negative N bias) (Sokolovskiy, 2003;Xie et al., 2006;Ao, 2007) which are shown as the grey curves in Fig. 1a and b.To correct the GNSS-RO refractivity retrieval bias in the presence of the ducting layer inside the PBL, a bilinear trapping layer model along with a reconstruction method were proposed by Xie06, which is described in the next subsection.

Bilinear trapping layer model and the reconstruction method
Xie06 demonstrated that an infinite number of refractivity solutions can generate the same bending angle observation using Eq. ( 2).Among all the refractivity solutions corresponding to the same bending angle profile, the one retrieved by the standard Abel inversion introduces the largest negative bias.
To mathematically describe each solution, Xie06 assume that the h(x) depicted in Fig. 1b can be approximated by a simple bilinear model, e.g., two connected straight line segments inside the trapping layer between h b and h t : Under this assumption, the missing information inside the trapping layer can be recovered by the parameterization of the h(x) between h b and h t , and the h 1 (x) at segment 1 can be derived analytically for given parameters x b , x m , h b , and h t : where h A (x) is the height function with respect to x from the standard Abel refractivity retrieval (grey curve in Fig. 1b) and  9), all the possible refractivity profile solutions, or the continuum of the solutions, below h t can be produced given five parameters defining the trapping layer: x b , x m , h b , h m , and h t .
To identify the best refractivity solution out of the continuum of the solutions, additional assumptions and constraints were proposed: A1.The height of the trapping layer top (h t ) is needed and can be derived from the standard Abel-retrieved refractivity profile.Due to the large bending angle (theoretically → ∞) near both the top and bottom of the trapping layer, the peak in the bending angle profile is assumed to be detectable and its corresponding impact parameter (x b ) can be identified accurately.The height of the trapping layer top h t , therefore, can be calculated in the Abel-inversion profile at the point with known x b .However, accurate detection of the trapping layer top is challenging in practice as it relies on the high-resolution bending angle profile which could be very noisy without filtering.On the other hand, the filtering process reduces the vertical resolution, which leads to error in the parameter x b .Therefore, a more robust method to detect x b parameter is needed.
A2.The standard Abel-retrieved refractivity profile near the top of the trapping layer behaves like a square-root function.Expanding by Taylor series around z = 0, the Abel-inversion result function h A (x) can be written as Keeping the leading order term in z, the function h A (x) − h t is assumed to behave like a parabola near x b and the derivative of its square can be written as so that C can be determined through linear regression in h A (x) over 200 m below x b (Xie et al., 2006).By neglecting higher-order terms of the derivative, C can be expressed as Since the parameters x b and h t are known, after C is determined the parameter x m can be calculated given h b using Eq. ( 12).
A3. Based on global observations of high-vertical-resolution sounding over the ocean, the slope of h 1 (x) is assumed to be continuous at h b to the bottom of the trapping layer: This assumption determines the parameter h m by linear extrapolation of h 1 (x) until x = x m .
Up to this point, all five parameters are connected so that once the parameter h b is given the other four can be determined.However, while a large number of h b values can be used in a single GNSS-RO case to generate a family of candidate profiles, choosing the correct profile from the family is still challenging.To determine the h b parameter one needs an additional constraint.
A4.The surface refractivity constraint is applied, which requires the extension of the RO bending angle observation to the Earth's surface, x 0 , i.e., from the whole family of refractivity profiles only the one with minimum height starting at zero will be chosen as the best solution of the reconstructed profile.However, for a number of reasons, including measurement noise and tracking error, horizontal variability, and diffraction effects, the retrieved profiles do not often reach the Earth's surface (Ao et al., 2012).Thus, the surface constraint is very difficult to fulfill and becomes the main obstacle in applying the reconstruction method in practice.
As shown above, the reconstruction method in Xie06 depends upon several conditions which may not be achievable.In this paper, we aim to refine this reconstruction method by improving the applicability of (A1) and replacing the surface refractivity constraint in (A4) with the combination of optimal estimation along with the external PW observation constraint.Furthermore, an improved method to determine the trapping layer top (x b ) is proposed and the bilinear parameterization model of the trapping layer is modified to reflect the smoother and more realistic structure of the h(x) curve.

Optimal estimation implementation
Our approach retains Eq. ( 9) from Xie06 as the core concept of the retrieval method.First, to improve the determination of the height of the ducting layer, a refined algorithm is developed, which no longer determines x b only by the noisy bending angle profile but also through optimal estimation iterations.Secondly, the retrieved PW from ancillary data is used to replace the surface constraint in (A4) (Xie06).This research focuses on correcting the N -bias below the trapping layer and deriving the key parameters such as duct altitude, thickness, and refractivity gradient information.
To implement these two major changes to the reconstruction algorithm, the basic parameterization process needs to be modified as described in the following subsection.

Parameterization
In the Xie06 reconstruction method, all five parameters are connected and only one free parameter, h b , is needed to define each refractivity profile.To relax the x b determination constraint in (A1), we use a different approach.
First, the highest 100 m of h 1 (x) will be replaced by the linear extrapolation from below; its slope is determined by linear regression between 100 and 200 m below the height h b .The reason is that the top of the analytical solution calculated by Eq. ( 9) is sensitive to the mis-modeling of h 2 (x) and h 3 (x), which are assumed to be bilinear segments between h b and h t .Spurious spikes and fluctuations are found at the top 100 m of the function h 1 (x) if the straight line assumption inside the critical layer is violated by the data or the parameters, e.g., h b , are not appropriately estimated from Equation (12).While these fluctuations in h 1 (x) are usually small (< 20 m) and typically occur within 100 m below h b , they can significantly change the slope at h b , which makes the h 2 (x) value obtained from assumption A3 inconsistent with having one refractivity for each impact parameter below h b .Therefore, we replaced the top portion of h 1 (x) with a linear extension of the curve below to remove the fluctuation.Note that the variable h b changes when the top of h 1 (x) is replaced.To avoid the inconsistency between the new h b and the originally specified h b , x m is chosen over h b as the "free" variable to construct a profile inside and below the critical layer.
Second, instead of assuming x b as known, both parameters x b and x m are used as "free" variables as is h b in the Xie06 approach.Namely, a pair of x b and x m are required to determine one refractivity profile in this new parameterization model.With both modifications, (A3) from Sect.2.1.2can be easily applied to find h b and h m by extending the top of h 1 (x) with the same slope until x = x b and x m , respectively (see Fig. 1b).These modifications, obviously, lessen the dependence upon (A1) but increase the reliance on the identification of multiple parameters x b and x m , which requires other constraints to assist in choosing the best candidate.In the next subsection we will investigate the use of PW observation as an additional constraint to choose the correct profile among the solutions.

Precipitable water as the external constraint
Precipitable water, PW, is the total column water vapor content in the Earth's atmosphere.In this research PW is chosen as the constraint over other physical quantities for several reasons.First, most of the water vapor in the atmosphere is located within the PBL so that accurate PW observations can provide extra information below the ducting layer to assist GNSS-RO retrievals.As an example shown in Fig. 2, each refractivity profile candidate corresponding to the same bending angle shows distinctive PW values.This is reasonable since the refractivity is strongly related to the water vapor content in Eq. ( 5) and the larger PW corresponds to greater refractivity.
Second, PW observations are available globally from a variety of sensors (Millan et al., 2016).The microwave radiometry remote sensing technology is used onboard several Earth-observing satellites including AMSR-E ( Kawanishi et al., 2003), AMSR-2 (Imaoka et al., 2010), TMI (Kummerow et al., 1998), and SSM/I (Alishouse et al., 1990) to acquire accurate PW observations over oceans.Also, the laserdiode sensor onboard airborne platforms such as AMDAR (Petersen et al., 2016) and ground-based GNSS receiver networks (Yuan et al., 2014) can provide PW observations over land.Here we focus on the AMSR-E observations onboard NASA's Aqua satellite, while PW observations from other sources can be applied to this method as well.AMSR-E conically scanned the Earth's surface with the water vapor band at ∼ 22 GHz and provided the precipitable water estimates over the oceans with limited error (∼ 0.6 mm) (Wentz and www.atmos-meas-tech.net/10/4761/2017/Atmos.Meas.Tech., 10, 4761-4776, 2017 Figure 2. The family of refractivity profiles calculated from a single simulated GNSS-RO bending angle profile (dotted line) using Eq. ( 9) and the parameterization method described in Sect.2.2.1.Each profile corresponds to a distinctive PW value, which can be used as a constraint for GNSS-RO retrievals within the boundary layer.Meissner, 2000).While the PW calculated by GNSS RO is negatively biased in the presence of ducting conditions, the collocated AMSR-E PW can be used to select the most optimal solution among the candidate profiles.Third, PW value of a candidate refractivity profile can be easily calculated to compare with the PW observations.In this study, the iterative direct method (Kursinski and Hajj, 2001) utilizing the temperature information from the European Centre for Medium-Range Weather Forecasts (ECMWF) atmospheric analysis is used to derive the PW from each candidate refractivity profile.The high-resolution ECMWF analysis data (TL799L91) used in this research have 91 vertical levels from the surface to 0.01 hPa and 0.25 • horizontal resolution.The data are modeled at every 6 h and unevenly sampled in vertical space which has higher resolution near the surface (∼ 40 m).The core concept of the direct method is derived by combining the hydrostatic and ideal gas laws.It was shown by Kursinski and Hajj (2001) that the relation between different atmospheric pressure levels can be described as where p is the pressure, T is the temperature, i is the index of each height interval, g is the mean gravitational acceleration, R is the universal gas constant, and z is height.m is the mean molecular mass of atmosphere which takes both dry air and vapor into account: where m v and m d are the molecular mass of dry air (∼ 28.97 g mol −1 ) and water vapor (∼ 18.02 g mol −1 ), respectively.Using Eq. ( 15) along with the refractivity Eq. ( 5) one can solve the water vapor pressure profile e iteratively by updating m at each step and the convergence at each height interval can be reached in one or two iterations.
The calculated water vapor profile, however, may not reach the surface of the Earth in most of the cases, while the PW calculation requires the information of water content in the atmosphere all the way to the surface.Hence, a reasonable assumption, which can be observed in many cases, is made that the moisture within the boundary layer is well mixed and the specific humidity q is constant from the surface up until the lowest point of the RO-retrieved profile.The specific humidity profile q(h) can be simply calculated (Stull, 2015) by By assuming the value of q at the surface equals the q(h) at the lowest height, we can calculate the precipitable water with the integration (Millan et al., 2016) from surface to the height when temperature reaches 230 K (Kursinski and Hajj, 2001): where the pressure profile p can be extrapolated from the lowest profile height to the surface using an exponential fit function.

Optimal estimation
To choose the best refractivity profile from a family of candidates with different x b and x m , an optimal estimation method (Rodgers, 2000) is used based on a Bayesian solution that minimizes the cost function of a linear inverse problem.In this method, the state vector s consists of two variables, that can be connected to the observation vector y with a forward model F where Because x b and x m are related, the second component of the state vector s is set as x m − x b instead of x m so that the two components can be treated independently.The observed PW is the observation vector y for the optimal estimation problem: Table 1.The spatial and temporal distance between different observation methods for each collocated case analyzed in this paper.The RAOB-AMSR column shows the differences between the RAOB and its closest AMSR-E measurement location in the simulation results, while the RAOB-RO column shows the differences between the actual RO tangent point location and its closest RAOB measurement.The RAOB observations in cases 1 to 3 are repeated in both the simulation and the actual data analysis.The forward model F to calculate the measurement y for each given state s is described in Sect.2.2.1 and 2.2.2.The Jacobian matrix K defined as is calculated numerically from the variation of y after perturbing the corresponding state s at the iteration step n.
Defining C s 0 and C y as the error covariance matrices of the a priori state s and the measurement y, one can estimate the best solution of s iteratively: where s 0 is the a priori guess of the state s and superscript T denotes the transpose of the matrix.For this study, the state a priori x b is determined by the impact parameter where a sharp transition occurs in the bending angle profile.The determination process using the step function correlation is described in Appendix A. However, the a priori information of the parameter x m − x b , which is highly correlated with the "strength" of ducting, cannot be obtained directly from the current GNSS-RO or AMSR-E measurement.In this study the x m − x b a priori value is chosen as the constant of 250 m, which is approximately the average number of x m − x b from all the radiosonde profiles (19 cases) used in this study.
To calculate the covariance matrix, the uncertainty of each variable is required.The uncertainty of x b is set as ±40 m, mentioned in Appendix A, and used to form the C s 0 matrix.The uncertainty of x m − x b , on the other hand, is not  1.
known and the a priori constant we chose is not based on any reliable sources.While the SD of x m − x b is ∼ 80 m in the radiosonde profiles, we conservatively set the uncertainty of x m − x b as large as ±400 m to allow the estimation of this parameter with more flexibility and insensitivity to the a priori constant we chose.The AMSR-E PW retrieval contains an error of ∼ 0.6 mm, but additional errors could rise from RO -AMSR-E collocation distances and forward modeling.Therefore, the conservative PW margin of 1 mm is used as the uncertainty of the PW observation in the C y matrix.Both C s 0 and C y are generated as simple diagonal matrices.Given appropriate initial conditions for all the least-squares fits included, the iterative process of Eq. ( 23) normally converges in a few iterations.The estimation results select the refractivity profile best fitted to the given x b , x m − x b a priori and PW observations, correct the N-bias, and provide the PBLtop information including its altitude and refractivity gradient.It should be emphasized that the optimal estimation also creates a framework for solving the ill-posed inversion problem of refractivity retrievals under ducting by incorporating multiple external constraints.In addition to PW constraint, the flexibility of the optimal estimation framework allows the use of other physical constraints to correct the N -bias in the presence of ducting.

Simulation results
To test and validate our algorithm, we conduct a simulation study utilizing radiosonde measurements from the VA-MOS Ocean-Cloud-Atmosphere-Land Study (VOCALS) campaign (Wood et al., 2011).The VOCALS campaign dataset is used because of its location in the southeastern Pacific Ocean, which has the world's most persistent subtropical stratocumulus deck at the top of the boundary layer (Bretherton et al., 2004).Moreover, the region also has one of the highest frequencies of ducting conditions (Lopez, 2009) and leads to a large N-bias in RO refractivity retrievals (Xie et al., 2010).Nineteen observations with strong refractivity gradient at the top of PBL are selected for the simulation, and six among them are collocated with AMSR-E measurements (Fig. 3 and Table 1).The x-h curves of these six cases are shown in Fig. 4. Original RAOB h(x) profiles are shown as the light red lines in Fig. 4, which are non-monotonic functions in the trapping layer near ∼ 1.5 km for all six cases.Using the RAOB refractivity profiles as reference, we generate an observed bending angle which is then Abel inverted to simulate the standard retrieved GNSS-RO refractivity profiles.While x is not monotonically increasing in the RAOB refractivity profiles, the forward calculation of Eq. ( 2) should be used here to generate the RO bending angle.Note that the potential errors caused by horizontal refractivity gradient are neglected in the bending angle simulation.The resulting standard Abel retrievals (x-h curves) are shown as black dotted lines.The Abel retrievals diverge from the RAOB profiles beneath the top of the ducting layer and cause negative bias in the x profiles below.The corresponding refractivity profiles for these six cases are shown in Fig. 5, where the standard Abel-retrieved RO refractivity profiles (dotted) contain large negative biases below the trapping layer when compared to the original RAOB profiles (light red lines).The collocated ECMWF analysis profiles are also shown as light green lines in Figs. 4 and 5.The ECMWF analysis tends to underestimate the ducting layer height and the refractivity gradient inside, which causes negative refractivity biases at lower altitude when compared to the radiosonde measurements.Since the VOCALS results were not assimilated in the ECMWF analysis, these two data sources can be regarded as indepen-dent.The statistically low PBL heights in ECMWF, which were extensively observed in the region, imply an erroneous refractivity profile below the ducting layer.This difference has been attributed to the model physics and assimilation process limitations (Xie et al., 2012).Even though ECMWF and other NWP (numerical weather prediction) systems assimilate both GNSS-RO bending angles and AMSR-E radiances, it is not clear that the full vertical resolution of the measurements can be taken into account.Thus, an independent, unbiased, refractivity retrieval outside of NWP data assimilation systems remains extremely valuable.
The biased GPS-RO refractivity simulations are then processed with the proposed optimal estimation method and compared to the original radiosonde refractivity profiles.The PW values calculated by Eq. ( 18) for each case from the RAOB profiles are shown in the lower-left corner of each panel along with collocated AMSR-E measurements.To validate our reconstruction method in this section, the background temperature and pressure profiles required in the direct method are given by radiosonde observations which are regarded as truth.The reconstruction results using the RAOB PW are presented in Figs. 4 and 5 as dark-red dashed lines.As Fig. 5 shows, when the unbiased RAOB PW values are used the reconstructed profiles can effectively reduce the Nbias and provide accurate estimates below the trapping layer.Optimal estimation also provides the estimated parameters x b and x m , which can be used to identify the altitude, thickness, and the refractivity gradient of the ducting layer.A discrepancy can be observed in the straight line section in the trapping layer corresponding to h 2 (x) and h 3 (x), where the original RAOB refractivity profile is not represented by two straight lines as we assumed.Fortunately, this approximation only induces small differences and has little impact on the reconstructed profile below the trapping layer.However, the h 3 (x) function inside the ducting layer is sensitive to the x b location and the slope at the top of the h 1 (x).This could lead to large refractivity differences from the true profile when the x b is not accurately determined or when the slope of h 1 (x) and h 2 (x) near x b are not continuous as expected.This error cannot be corrected without adding additional constraints or measurements to further determine the vertical structure inside the trapping layer.
In practice, the PW information from RAOB is not always available for nearby GNSS-RO soundings due to the sparsity of radiosonde stations in remote areas.To demonstrate the ability of using other ancillary PW sources in the proposed algorithm, the reconstructed profiles using the collocated ECMWF and AMSR-E PW are also presented in Fig. 5 as dark-green dashed lines and blue dashed lines, respectively.All three different PW values used in the reconstruction method are listed in the lower-left corner of each panel.
The PW values acquired from the three external sources (RAOB, ECMWF, AMSR-E) in all six cases are greater than the ones calculated from the negatively biased Abel-inverted profiles, which suggests dry biases in the Abel retrievals inside the boundary layer when ducting occurs.Therefore, the reconstructed profiles from the optimal estimation with larger external PW should lead to larger refractivity inside the PBL and mitigate the N -bias.
The statistical results of the 19 RAOB cases using the reconstruction method with the radiosonde PW are shown in Fig. 6a.The refractivity difference is defined as where N RAOB is the radiosonde refractivity and N RO is the standard Abel refractivity retrievals (dotted lines) or the reconstructed profiles (red lines) in Fig. 6a.As depicted in Ao (2007), the negative N-bias reaches the greatest value (−8 to −17 %) near the height of h t and decreases to ∼ −5 % near the Earth's surface.On the other hand, using the optimal estimation method the retrieved refractivity profiles remain unbiased (< 1 %) below the trapping layer.The large error up to +13 % at the top of the ducting layer (h 3 (x)) is due to the presence of the sharp refractivity gradient, and minor x b and h m estimation differences from the RAOB profiles could lead to large differences.While the error mostly occurs around the Figure 6.The refractivity differences between the RAOB profiles, the Abel-retrieved profiles, and the reconstructed profiles for the 19 simulation cases with the single ducting layer in the VOCALS campaign using the PW from RAOB (a) and ECMWF (b).While the proposed method can correct the N -bias using both RAOB and ECMWF PW, the reconstructed profiles show higher variance and are slightly biased below the trapping layer with ECMWF PW.
top of the ducting layer, it has very limited effects on the estimated profile below and the character inside the trapping layer.However, the errors in external PW constraints will affect the reconstruction results.As presented in Fig. 6b, while the reconstructed results using ECMWF PW reduce the N-bias, it still leads to a smaller negative bias in the reconstructed results (−1.54 % on average) compared to the ones reconstructed from the RAOB PW (−0.01 % on average).This may be due to a systematic underestimation of PW by the ECMWF analysis.Approximately 1 mm of the PW bias can cause a ∼ 3 % refractivity bias at the height h t and ∼ 1 % at the surface.Although the slight negative bias caused by lower PW values (∼ 1 mm) could reduce its reliability, these results suggest that the ECMWF analysis can still be used to improve the retrieval under the trapping layer.
The statistical results of the six cases using the three different PW sources are shown in Fig. 7.While the reconstruction retrievals using ECMWF PW are negatively biased below the trapping layer, the results using the collocated AMSR-E measurements tend to be non-biased in general because of relatively better agreement between AMSR-E and radiosonde PW.A positively biased outlier in AMSR-E PW reconstruction results can be identified as case 5, whose AMSR-E PW is apparently larger than those from other sources.The AMSR-E, ECMWF, and RAOB PW value comparisons are presented with the scatter plot in Fig. 8, where ECMWF PW shows a clear negative bias compared to the RAOB PW.The large difference of PW in case 5 and case 1 may be due to the large distance (431 km) and the temporal difference (1.5 h) between the AMSR-E measurement and the RAOB location, respectively.However, the cause of the PW difference in case 4 is unknown.While more analysis is needed to identify the true cause of the limited bias in each individual case, the microwave PW retrieval is still proven to be another useful constraint when closely collocated with GNSS RO and provides a feasible solution for GNSS RO during ducting.

Actual GNSS-RO data results
We now apply our reconstruction method on actual Constellation Observing System for Meteorology, Ionosphere, and Climate (COSMIC) RO data.Eight COSMIC occultations collocated with VOCALS radiosondes (Fig. 9) are chosen.Three criteria are utilized for choosing these cases: a spatial distance of less than 300 km and a temporal difference of less than 3 h -the lowest height of the GPS-RO refractivity profile reaches below 1 km to ensure the trapping layer is included.We also exclude the cases with complex x-h structure inside the trapping layer, which can heavily violate the bilinear assumption, and the cases with multiple ducting layers, which makes Eq. ( 9) inapplicable.Approximately 15 % of the total number of cases are ruled out by these two additional requirements.The first three GPS-RO cases have both collocated radiosondes and AMSR-E measurements and are numbered as cases 1 to 3 in Table 1.The other five cases which do not share the collocated RAOB measurements are numbered 7 to 11 for clarity.In practice, RAOB temperature and pressure profiles at the GNSS-RO collocation may not be available for the RO's PW calculation when using the direct method described in Sect.2.2.2.Therefore, in this section we also include the ECMWF analysis profiles to compute PW for the GNSS RO and compare them with corresponding radiosonde and AMSR-E measurements.
The results are shown as dashed lines in Fig. 10.Similar to the simulation results in the previous section, the actual COS-MIC RO refractivity profiles (dotted lines) are negatively biased compared to the collocated RAOB and ECMWF analysis.This negative refractivity bias leads to a smaller GPS-RO PW value than the one calculated from the given RAOB profiles, ECMWF profiles, and AMSR-E measurements.The reconstructed results correct the bias to different degrees based on the source of PW ancillary data.Two main differences can be observed when comparing the reconstructed profiles to the reference collocated radiosonde profiles.First, the refractivity profile above h t from RO and radiosondes is not exactly the same, in contrast to the assumption that the bias only comes from ducting.This can bias the PW value calculation for all possible candidate profiles.Second, the estimated x b in reconstruction results can have at most a 200 m difference with the corresponding radiosonde observations and cause the different shape of the reconstruction even if the PW is obtained accurately.These two differences can be caused by the spatial-temporal difference between RO and RAOB observations, which are normally more than 200 km and 1 h apart.For example, in the reconstructed profile for case 3, COSMIC sounding is only 10.7 km and 1.18 h apart from the RAOB location, which agrees well with the RAOB profile in x b and the refractivity profile above.Another possible cause of x b discrepancy is the error in GNSS-RO measurement due to horizontal inhomogeneity in the atmosphere and the ionosphere (Zeng et al., 2016).In ducting conditions, this error can be amplified and can shift the impact parameter of the top of the boundary layer for more than 100 m.While addressing the horizontal inhomogeneity is beyond the scope of this article, the impact of the horizontal refractivity gradient on the reconstruction method can be further investigated in future work.
The statistical results of the refractivity difference compared to the collocated radiosondes are shown in Fig. 11.Two main differences mentioned in the previous paragraph can be easily seen: the differences above the ducting level can be as much as 5 %, and the x b difference cause < −10 % error above h t and > 10 % error under h t .The not-so-closely Figure 10.Refractivity profiles from the actual RO (dotted lines), collocated radiosonde measurements (solid red lines), ECMWF analysis (solid green lines), and reconstruction using RAOB (red dashed lines) and ECMWF-computed PW (green dashed lines) for the eight collocated cases in the VOCALS campaign.The reconstructed profiles using collocated AMSR-E measurements are also shown (blue dashed lines) in the first three cases.All case numbers are in the left side of each panel.Although the ducting layer heights in most cases are different from the RAOB profiles due to the distance and different footprints between the two observations, the RO reconstruction results are still able to correct the N-bias below the trapping layer with a higher amount of water vapor content measured by the AMSR-E or ECMWF model.
Figure 11.The refractivity differences of the Abel-retrieved and the reconstructed profiles from different PW sources compared with the original RAOB profiles in the eight actual data cases.Compared to the 15 % negative bias below the height h t from the Abel-inversion results, the reconstructed profiles utilizing closed PW sources can limit the error to within 5 %.The large negative bias and variance above h t are mostly due to the spatial and temporal distances between the RO and RAOB, which are normally more than 200 km and 1 h apart.As the results shown in the simulation, negative biases can still be observed in the ECMWF PW reconstructed profiles in actual data cases.collocated RAOB profile PW (red lines) can still maintain the N -bias below the trapping layer at less than 5 % in all the cases without bias.The three cases using AMSR-E PW retrievals shown in blue lines agree better with the reference RAOB profiles, which can be attributed to the unbiased AMSR-E PW observation.Like the simulation results, the reconstructed profiles using ECMWF PW are negatively biased (∼ −2 %) against the ones using other PW sources.This is because the PW calculated by ECMWF is negatively biased (Fig. 8).Overall, our results show that reconstructed profiles utilizing external PW sources can substantially reduce the negative N -biases and limit the error to within 5 % with zero mean below h t from the 15 % negative error in the standard Abel refractivity retrievals.

Conclusions
GNSS RO has been extensively used in atmospheric profiling and weather forecasting.However, the RO profiling in the presence of the ducting layer remains a challenge.Ducting, which occurs when the refractivity gradient exceeds the critical refraction (−157 N-unit km −1 ), can cause the lack of bending angle information within the trapping layer.As the retrieved bending angle loses its one-to-one relationship with the atmospheric refractivity, the standard Abel inversion gives the refractivity solution with the largest negative bias.By approximating an analytical solution to the profile below the trapping layer and introducing a series of constraints, the method by Xie06 is able to reconstruct the pro-file based on GNSS-RO observations.However, the reconstruction method in Xie06 relies on several idealizations that are difficult to implement considering the uncertainty of the real RO measurements.To develop a practical reconstruction method, this paper validated a new implementation framework to incorporate constrained optimal estimation into the RO retrievals in the presence of a ducting layer.
The proposed method modified the parameterization process to include more free parameters and reduce the reliance on idealized assumptions.The optimal estimation method is used to select the candidate that minimizes the cost function, which is defined by the difference between the known reference (i.e., ancillary PW observations and a priori state) and those calculated from each retrieval.PW observations, which can be obtained by remote sensing instruments such as AMSR-E, can serve as an external constraint in the reconstruction method.The process to infer the boundary layer height from bending angle profiles has also been refined to provide a robust and accurate estimation of a priori x b .The new reconstruction method has been applied to both the simulated GNSS-RO profiles and actual GNSS-RO data.The results show that given accurate PW the proposed method greatly reduces the reconstruction error to less than 1 % in simulation and 5 % in actual cases.While the method cannot fully reconstruct the vertical structure inside the trapping layer, the iterated parameters are able to give improved estimation of PBL-top features including ducting layer altitude, thickness, and the refractivity gradient.
To improve this reconstruction technique, several sources of uncertainty need to be further examined.The biases in different PW sources should be identified before being used as constraints, and the impact of the spatial and temporal difference between the chosen PW observations and GNSS RO requires further investigation.Also, the deviation from the assumptions of constant specific humidity from surface up to the minimum height of RO sounding and the continuity of x(h) slope at the bottom of trapping layer may cause additional errors in the reconstructed profile.The optimal estimation method developed in this work can be improved by incorporating other potential observations and constraints in the future and will help to better characterize the vertical structure of the PBL globally using GNSS-RO measurements.It should be recognized that the absolute accuracy of the reconstructed GNSS-RO refractivity is influenced by the uncertainty of the external constraints.The lower SI traceability of the reconstructed refractivity within the PBL compared to the upper troposphere and lower stratosphere (UTLS) region can limit its applicability in long-term climate monitoring.

Appendix A: Determination of x b in the bending angle profiles
In this appendix we describe a new method to detect the impact parameter x b where the ducting occurs.Theoretically, the bending angle should reach infinity when the tangent point of the signal path is located inside the trapping layer (Sokolovskiy, 2003).In practice, the infinite value of bending angle is not observable in a finite observation, but the singularity in the bending angle close to the impact parameter x b results in a sharp transition.The location of the sharp transition of the bending angle measurement can provide us with valuable information on the a priori x b .In this paper, the parameter x b is determined by the peak of the correlation between the high-resolution bending angle profile and a step function using a two-step approach, which is similar to the wavelet covariance transform method proposed by Ratnam et al. (2010).The step function we used has the value of +1 at its lower 500 m and −1 at its higher 500 m, which is similar to the shape of the bending angle profile affected by ducting.The altitude of the transition from −1 to +1 in the step function matches the sharp transition of the bending angle at the ducting layer.The pattern matching correlation result is shown as a blue solid line in Fig. A1a.In this figure, the highresolution bending angle is shown in grey and shows a sharp transition around the impact parameter a = 6367 km.Note that the 1 m resolution bending angle is used instead of the common low-resolution profile which has been filtered with a 200 m window and degraded x b precision.The correlation shows a clear peak because the 1 km length step function filters out most of the fluctuations caused by noise, multipath, or highly variable water vapor content close to the ducting layer.The maximum of the correlation function indicated by the dashed line is close to the impact parameter where the sharp transition of the bending angle occurs.
While this peak can provide the coarse estimation of x b within 250 m, the length of the step function is very insensitive to the transient behavior of the bending angle.To enhance the precision of the x b estimate, the second correlation with a shorter step function is used.In this search, a 150 m length step function is used to repeat the correlation with the high-resolution bending angle profile.However, when the tangent point lies close to the top of the ducting layer, the determination of the sharp transition becomes difficult due to fluctuations in the bending angle.On the other hand, the observed bending angle rarely contains large fluctuations below the trapping layer.Therefore, to put more weight in correlation at the lower half of the step function, it has been modified into an asymmetric shape, with a value +1 for the lower 90 m and a value −1 for the higher 60 m, which extends the lower part while shrinking the upper part.In addition, the bending angle has been de-trended from the exponential fitting function before applying the second correlation to simplify the profile which can focus on the transition due to ducting instead of normal refractivity increases.The result of the sec- ond correlation is shown in Fig. A1b as the green line.Due to the shorter integration period the second correlation has a higher variability than the first one.The peak of the second correlation is then searched for over the range between −250 and 250 m relative to the maximum of the first correlation function that covers the impact parameter where the sharp transition could occur.As Fig. A1b shows, the second correlation peak location, shown as a green line, can clearly determine the location of the sharp transition to be used as the estimated x b .Using this method, x b can be determined with an uncertainty of less than 50 m.

Figure 1 .
Figure 1.The illustration of the corresponding refractivity profile (a) and the h(x) function (b) when ducting occurs.The true profiles are shown in black and the profiles acquired by GNSS-RO Abel retrievals are shown in grey.The large refractivity gradient between h m and h t causes the negative slope of h 3 (x) and the multivalued function h(x) between x b and x m .

Figure 3 .
Figure 3.The map of the six collocated RAOB and AMSR-E measurements in the VOCALS campaign.The spatial and temporal differences for all six cases are listed in theTable 1.

Figure 4 .
Figure 4.The x-h relationship of simulated RO, radiosonde measurements, ECMWF analysis, and reconstruction using RAOB and ECMWF-computed PW for the six collocated cases.The case numbers are shown in the lower-right corner of each panel.As shown in the figures the RO simulations maintain a one-to-one relationship between x and height when ducting happens, which causes negative bias compared to the RAOB results.The proposed method can reconstruct the bilinear shape inside the trapping layer and correct the N-bias below.

Figure 5 .
Figure5.Refractivity profiles from simulated RO, radiosonde measurements, ECMWF analysis, and reconstruction using RAOB and ECMWF-computed PW and AMSR-E retrievals for the six collocated cases.The case number is found on the left side of each panel.It can be observed that the ducting layer height in the ECMWF model is mostly lower than the one measured by radiosondes.The RO optimal estimation results, which correspond to different PW sources, can correct the N -bias with a higher amount of water vapor content measured by other techniques.

Figure 7 .
Figure 7. Refractivity differences of the Abel-retrieved and the reconstructed profiles from different PW sources compared to the original RAOB profiles in the six simulated cases.It can be seen that the results using ECMWF PW are negatively biased (−1 to −5 %) while the collocated AMSR-E reconstructed results more closely aligned with the RAOB profiles.An outlier of the AMSR-E reconstruction shown in the figure is case 5, which was measured 431 km away from the corresponding RAOB case.

Figure 8 .
Figure 8. PW scatter plot for simulation results from different sources.All 19 RAOB cases PW are compared with the collocated ECMWF and AMSR-E PW in this figure.The x axis is the RAOB PW, and the black line is the identity line.ECMWF PW values are systematically lower (1 to 2 mm) than RAOB PW and cause negative biases in reconstructed profiles.

Figure 9 .
Figure 9.The map of the eight collocated RAOB and COSMIC GPS-RO measurements in the VOCALS campaign.Three AMSR-E measurement locations, which coincided with RAOB measurements in the first three cases, are shown in blue squares.The distance between the RAOB (and AMSR-E) and the corresponding RO cases is 250 km on average.The temporal difference for all eight cases is within 3 h.

Figure A1 .
Figure A1.Ducting layer height determination using bending angle profiles.The correlation result with the long step function is shown in (a) with a blue line, which has been scaled and shifted for demonstration.The peak location identifies the approximated ducting layer height.The correlation result with the short step function is shown in (b) with a green line.Panel (b) is the enlarged image of the dashed-line rectangle shown in panel (a).The corresponding impact parameter of the sharp bending angle transition can be found at the location of the second correlation peak within the range of the first correlation hump (±250 m).The correlation results in both panels are scaled down and shifted to fit into the figures.