GNSS-RO Refractivity Bias Correction Under Ducting Layer Using Surface-Reflection Signal

Abstract

Due to its high vertical resolution and cloud-penetrating capability, GNSS-Radio Occultation (RO) remote sensing technique has been utilized to observe the vertical structure of the Planetary Boundary Layer (PBL) in recent years. However, the critical refraction, or ducting, caused by large refractivity gradients usually associated with the top of the stratocumulus clouds, can negatively bias the retrieved refractivity and humidity within the PBL. Previous research has shown that combining RO retrievals and the external information, such as collocated precipitable water (PW) estimates, can effectively reduce the negative bias and enhance the retrieval quality. Nevertheless, the requirement of collocated observations from other techniques limits the applicability of this reconstruction method in practice. Here, we describe an alternative approach that uses the coherent grazing signals from the same RO event that are reflected by the Earth’s surface to remove the bias due to ducting. Additional observations are no longer necessary in this approach because the reflected signals provide the extra constraint. A least squares framework is used to select the candidate from a family of solutions wherein reflected bending angle best matches the corresponding observation. This new method was validated by both multiple phase screen (MPS) simulation and the simulated RO bending angle via forward Abel transform, and it was tested with the actual GPS-RO measurements. While, in general, the reflected bending retrieved from the current mission was noisy, the results show that this approach can potentially reduce the negative bias and improve RO observation within the PBL without assistance by the external information, such as PW.

1. Introduction

Global Navigation Satellite System Radio Occultation (GNSS-RO) is an atmospheric limb sounding technique to observe the vertical thermodynamic structure using GNSS signals received by low Earth orbit (LEO) satellites. When the occulted signal from the GNSS constellations propagates through the stratified atmosphere, its ray path could bend at the limb of the Earth due to the change of refractivity [1]. The bending angle of each ray path, which can be calculated by the time-varying signal phase delay observations, will be inverted to retrieve the vertical profiles of refractivity, temperature, and moisture. Since the first GNSS-RO mission, Global Positioning System/Meteorology (GPS/MET), successfully demonstrated the remote sensing technique in 1995, several new GNSS-RO instruments, including CHAMP, SAC-C, COSMIC, GRACE, and METOP, were launched and currently provide thousands of daily profiles in combination. Recently, launched GNSS-RO missions, such as COSMIC II (June 2019) [2], METOP-C (November 2018), and GRACE-FO (May 2018), were also deployed to increase the total number of RO atmospheric observations to more than 10000 per day and further improve the global weather prediction.

In recent years, GNSS-RO remote sensing has been applied to planetary boundary layer (PBL) studies [3,4]. PBL, which is the bottom layer of the atmosphere (≲ 2 km), controls the energy distribution from the surface and has a great impact to both local weather and global climate through turbulent winds and cloud formation [5]. Because of its importance in weather modeling, PBL has been extensively studied for decades and is listed as one of the most targeted observables in the decadal survey [6]. GNSS-RO technique is valuable to PBL probing tasks for several reasons. First, as a spaceborne remote sensing technique, GNSS-RO can provide worldwide boundary layer information, including remote regions, with few in-situ measurements. Second, GNSS-RO provides high vertical resolution (∼100 m) [7] refractivity measurements, which are not available from other nadir remote sensing instruments [8]. Third, the L-band GNSS signals, which are designed for navigation purposes, can penetrate clouds and precipitations [9], which are common at the level of PBL height. These advantages make GNSS-RO a unique and valuable technique for PBL characterization [10,11,12,13,14].

However, the strong temperature inversion and moisture change at the top of PBL could induce large negative refractivity gradient and sharply increase the bending of the propagating radio signal ray path. When the bending angle is larger than the curvature of the Earth’s surface, ducting occurs and the signal will be "trapped" inside the atmosphere. In general, ducting occurs when dN/dr ≲−157 (N-units/km) (the region between $h_m$ to $h_t$, as shown in Figure 1), which can be frequently observed in the subtropics below 2 km. Contrary to the ducting layer definition from von Engeln and Teixeira [15] ($h_m$ to $h_t$), here, we define it as the layer between $h_b$ and $h_t$, where the radio signals within will theoretically be trapped [16]. In the case of GNSS-RO, in which both the transmitter and receiver are located outside the atmosphere, the signal path will never be "trapped" inside the ducting layer. Instead, the tangent points of RO ray paths skip the ducting altitudes and cause the bending angle information loss inside the ducting layer. As a result, the refractivity retrieval process becomes an ill-posed inversion problem (i.e., multi-refractivity solutions correspond to the same bending angle profile) when ducting is present. The refractivity profile retrieved by standard Abel inversion, which assumes no ducting occurs and no bending angle gap within the ducting layer, will be negatively biased (N-bias) below the ducting ($h_t$), as shown in Figure 1 [16,17,18]. Note that only the elevated ducts ($h_b>0$) condition could cause the N-bias, whereas in the surface and evaporation duct cases, the RO profile would be cut-off at $h_t$, resulting in higher penetration altitude. Ao [16] and Xie et al. [18] stated that the N-bias caused by the ducting can be as much as 5% in several subtropical ocean basins. Therefore, it is essential to remove the N-bias in the retrieval profiles due to ducting for correctly characterizing PBL using GNSS-RO.

It is worth clarifying here that ducting is not the only cause of the GNSS-RO N-bias. Tracking error [19], cycle slips [20], unbalanced noise spectrum [21], and atmospheric turbulence [22] could also lead to both bending angle and N-bias in the lower-altitude RO retrievals. It is difficult to identify and distinguish different sources of N-bias solely by RO observations, and in this study, we only focus on correcting the N-bias caused by ducting. Several efforts were made to reduce the negative bias. Xie et al. [18] analytically proved that, when ducting occurs, one bending angle profile w.r.t. impact parameter can correspond to an infinite number of refractivity profile solutions. These refractivity solutions can be described by a single non-linear function which depends on the thermodynamic structure within the ducting layer. While the structure inside the ducting layer cannot be directly known from the GNSS-RO signal observations, one can use physical constraints to identify the correct profile from the solution continuum. The approach proposed by Xie et al. [18] constrains the solution by assuming GNSS-RO profiles can reach all the way down to the Earth’s surface and cut-off right at the surface level. In practice, however, the lowest point of GNSS-RO bending angle profiles do not often reach the surface due to measurement noise and tracking error [10], which makes the application of this constraint difficult, in practice. As an alternative constraint, Wang et al. [23] proposed distinguishing different profiles in the continuum with the integrated precipitable water (PW) and using the collocated external PW retrievals from other instruments as the constraint. This approach requires collocated observations, which are not always available and could induce the collocation error in the N-bias correction process.

In this paper, we introduce the concept of using reflected GNSS-RO signal for N-bias removal. Reflection component of the GNSS-RO signal can be easily detected (more than 45% of the cases) in the ocean basins with latitude greater than 20 degrees [24]. While the bending angle of the direct signal for each N profile in the continuum is the same, the corresponding reflected bending angle is distinctive and can be used to discriminate different profiles when ducting occurs. This concept will be further explored in Section 2 with Multiple Phase Screen (MPS) simulation, along with the reflected bending angle calculation. Since the reflection component can be extracted directly from the standard RO observations, no external information is required to set up a constraint. In Section 3, the least square approach using reflected bending angle as the constraint for choosing the unbiased refractivity solution is illustrated. The retrieval results in end-to-end simulation and actual cases will be shown in Section 4. Some restrictions of applying reflected bending as the constraint are discussed in Section 5. The conclusion is provided in Section 6. All the variables used in the main text are listed and defined in

Table 1. Symbols used throughout the main text of the article.

Symbol	Definition	Symbol	Definition
${\alpha }_D$	Direct bending	${\alpha }_R$	Reflected bending
${\alpha }_R^C$	reflected bending candidate	${\alpha }_R^O$	Observed reflected bending
${\Phi }_D$	Direct signal phase	${\Phi }_R$	Reflected signal phase
$\varphi$	grazing angle	${\rho }_s$	Reduction factor
a	Impact parameter	$a_S$	Impact parameter at the surface
$A_D$	Direct signal amplitude	$A_R$	Reflected signal amplitude
$d_{th}$	Ducting layer thickness	g	Rayleigh roughness parameter
h	Height	$h_b$	Height of ducting layer (bottom)
$h_m$	Height of ducting layer (middle)	$h_t$	Height of ducting layer (top)
$h_{cos}$	Cosine window width	i	Level index number
N	True refractivity	$N_{AI}$	Abel-inverse refractivity
$N_{RAOB}$	Radiosonde refractivity obs.	$S_h$	Surface standard deviation
n	True refractive index	$n_{AI}$	Abel-inverse refractive index
r	Radius	$r_E$	Earth’s radius
s	Total signal	$s^{\prime }$	Frequency shifted signal
t	Time	w	Window function
x	$n\left(r\right)r$	$x_b$	x value at bottom of ducting layer
$x_m$	Local maximum of x	$\lambda$	Signal wavelength

.

2. Reflecting Signal and Ducting

To illustrate that refractivity profiles are under-determined from the bending angle measurements in the presence of ducting, a simple example using Multiple Phase Screen (MPS) simulation is shown in Figure 2. A ducting layer located at 2 km can be simulated using a nonlinear function [25], as the blue line is shown in Figure 2a:

$$N\left(h\right)=350exp\left(\frac{-h}{7}\right)\times \left[1-\frac{0.2}{\pi }arctan\left[\frac{h-2}{0.06}\right]\right],$$

(1)

where N is the refractivity (N-unit) as a function of the height h (km). The minimum gradient $dN/dh$ can reach −250 N-units/km, which is less than the critical ducting value −157 N-units/km. The corresponding bending angle $\alpha$, with respect to the impact parameter a from the direct signals, can be calculated using a generalized forward Abel transform [17,18] to integrate the multi-value function $n\left(a\right)$:

$${\alpha }_D\left(a\right)=-2a\left[{\int }_a^{x_m}+{\int }_{x_m}^{x_b}+{\int }_{x_b}^{\infty }\right]\frac{1}{n\left(x\right)}\frac{dn\left(x\right)}{dx}\frac{dx}{\sqrt{x^2-a^2}},$$

(2)

where $x_m$ is the local maximum of x, and $x_b$ is the critical value when ducting starts occurring. The subscript D is added in Equation (2) to identify it as the bending angle of the direct ray path. The forward transformed bending angle profile, shown in Figure 2b, simulates the RO bending angle observation derived from direct GNSS signal paths. Note that the bending angle increases sharply at the height where the ducting occurs. One can apply the standard Abel-inversion operation to the simulated bending angle:

$$n_{AI}\left(x\right)=exp\left[\frac{1}{\pi }{\int }_x^{\infty }\frac{{\alpha }_D\left(a\right)da}{\sqrt{a^2-x^2}}\right],$$

(3)

$$N_{AI}(h)=[n_{AI}(h)-1]\times {10}^6,$$

(4)

where $n_{AI}$ is the Abel-inverse refractivity. Since the Abel-inversion assumes the variable $x\left(r\right)=nr$ is a single-valued function and the multi-valued portion in the profile is ignored, the resulting refractivity $N_{AI}\left(h\right)$ could contain a large negative bias compared to the true profile $N\left(h\right)$, as shown in Figure 2a. The fact that both refractivity profiles correspond to the same bending angle indicates that using only the direct GNSS RO signal cannot distinguish them in the retrieval process when ducting occurs. Unlike the approach that relies on the external information provided by other techniques to distinguish an individual profile, the feasibility of using a GNSS-RO signal itself that reflected from the the Earth’s surface was assessed in this research.

Here, we simulate the reflected GNSS-RO signals by following the steps of Beyerle et al. [26] to perform the MPS simulation with the refractivity profiles shown above. MPS modeled the atmospheric refractivity as a series of screens (as shown in Figure 3), and the propagating GNSS signal amplitude and phase through the screens can be calculated by iterating Fourier Transform and inverse Fourier Transform at each screen [25]. In our approach, there are $2^{19}$ grid points for each screen that expands from $-274.288$ km to 250 km vertically, as well as 4001 phase screens in total from $-1000$ km to 3000 km horizontally. The interval between each grid is 1 m, and the horizontal distance between each screen is 1 km. The radius of the Earth is set as 6370 km. In the last (observation) screen, we assume the vertical speed of the satellite is $2.7$ km/s corresponding to the satellite with orbit altitude $\sim 450$ km [26].

The resulting complex signal $s\left(t\right)$, which is the combination of direct and reflected signals, can be written as a function of time t:

$$s\left(t\right)=A_D\left(t\right)exp[i{\Phi }_D\left(t\right)]+A_R\left(t\right)exp[i{\Phi }_R\left(t\right)],$$

(5)

where $A_D$ and $A_R$ are the direct and reflected sub-signal amplitude, and ${\Phi }_D$ and ${\Phi }_R$ are the direct and reflected sub-signal phase, respectively. By adding a cosine window function above the Earth’s surface for each screen in the MPS simulation:

$$w\left(h\right)=\left\{\begin{array}{cc}1\hfill & :h>h_{cos}.\hfill \\ co{s}^2\left(\frac{\pi }{2}\frac{h_{cos}-h}{h_{cos}}\right)\hfill & :0<h<h_{cos}.\hfill \\ 0\hfill & :h<0.,\hfill \end{array}\right.$$

(6)

one can remove the reflected signal component (second term in Equation (5)) from s(t) and extract ${\Phi }_D\left(t\right)$. The variable $h_{cos}$, which is the layer width where the cosine window applies, is set as 2 km in this research. With this direct signal phase information, the spectrum of the signal s(t) can be shifted by the following equation and centered on the direct frequency, which is shown in Figure 4.

$$s^{\prime }\left(t\right)=\frac{s\left(t\right)}{exp[i{\Phi }_D\left(t\right)]}=A_D\left(t\right)+A_R\left(t\right)exp[i({\Phi }_R\left(t\right)-{\Phi }_D\left(t\right))].$$

(7)

As shown in Figure 4a, the signal consists of several frequency components due to different ray paths in the presence of ducting. The component located at zero-frequency goes through the direct ray path (Path 1 in Figure 4b), and the lower amplitude component shown before 85 s is the reflected sub-signal (Path 4 in Figure 4b). Because of the 50 Hz sampling rate, the reflected sub-signal is aliased before 50 s and gradually converges with the direct signal from 50 to 85 s, as a result of the closer ray path length between direct and reflected signals in the low elevation geometry. The direct signal vanishes at about 80 s because of strong defocusing at the ducting layer, but the signal component from the path below the ducting layer (Path 2 and 3 in Figure 4b) takes over in a relatively short time (∼82 s). This component, which starts at the frequency of ∼−2Hz, breaks into two different branches due to the strong refractivity gradient of ducting. One is closer to the ducting layer (Path 2 in Figure 4b) and dominates the received RO signal (located at zero-frequency) from 85 s to the end of the occultation period. The other one approaches to the Earth’s surface (Path 3 in Figure 4b) and eventually combines with the reflected component at around 85 s, when the tangent points of both sub-signal ray paths reach the surface and stop tracking. It should be emphasized that the tangent points of the reflected signals are always located below the surface, which make the reflected signals an effective constraint since no information gap will be formed due to ducting.

In this paper, the reflected GNSS-RO signals are utilized by calculating their bending angles. The bending angle of a reflected RO signal path (subscript R) based on a known ducted refractivity profile can be defined as [24]:

$${\alpha }_R\left(a\right)=-2a\left[{\int }_{a_S}^{x_m}+{\int }_{x_m}^{x_b}+{\int }_{x_b}^{\infty }\right]\frac{1}{n\left(x\right)}\frac{dn\left(x\right)}{dx}\frac{dx}{\sqrt{x^2-a^2}}-2arccos\left(\frac{a}{a_S}\right),$$

(8)

where $a_S$ is the impact parameter at the surface. The first term of Equation (8) describes the ray path bending due to the ducted atmosphere above the surface, while the second term is caused by the reflection of the surface. Note that the integration interval starts from the surface instead of the regular ray path impact parameter, which in the case of reflection, is located below the surface. Here, we focus on the reflections from a horizontal surface (e.g., ocean surface) without tilting so that the reflection points shifting described in Beyerle et al. [26] is not considered in this analysis. The roughness of the sea surface has limited impact in this case since the grazing angle is relatively small ($\sim 1^{\circ }$ when the impact height is 1 km below the surface).

The bending angle of each reflected path can be directly calculated from the received signal amplitude and phase via the current radio-holographic method, like phase matching (PM), by extending the impact parameter range below the surface [27]. Applying PM to the MPS, simulated signals result in both direct and reflected bending angles for the corresponding refractivity profile. In this analysis, the two refractivity profiles described by Equation (1), (4) are used in separate MPS simulations, and both direct and reflected bending angle for each refractivity profile are calculated by PM. The results are shown in Figure 5. The portion of the bending angle where the impact parameter range is above the surface (non-shaded) is contributed by the direct GNSS-RO signals, while the part below the surface (shaded) comes from the reflected GNSS-RO signal components only. Two bending angle profiles corresponding to the refractivity profile of Equation (1) (blue line) and Equation (4) (orange line) align to each other above the surface, which is consistent with our former analysis, shown in Figure 2. However, a large difference between these two reflected bending angle profiles can be found below the surface. The difference can reach 0.005 rad in reflected bending or 0.3 km in the impact parameter close to the surface and gradually decreases below. This difference in the reflected bending angle profile can thus be used to distinguish different profiles corresponding to the same direct bending angle in the ducting situations. Note that this is a useful feature for weather forecast data assimilation, which can identify or flag a specific case if it is contaminated by the ducting-caused N-bias [28]. Here, we take this idea further to reconstruct the refractivity profile and reduce the N-bias, and the algorithm to identify the correct profile with the reflected bending angle is introduced in the next section.

3. Bias Correction

3.1. Family Solution and Reflected Bending Angle

Xie et al. [18] proposed that, when ducting occurs, an infinite number of refractivity profiles can correspond to a single bending angle observation which caused an ill-posed inversion problem. Given an Abel-inverse refractivity profile, one can analytically derive the “family" of the candidate profiles $N_f\left(h,dx\right)$ [16,18], as shown in Appendix A, where $dx=x_m-x_b$. While, in Wang et al. [23], two independent variables were used to account for the uncertainty of $x_b$, here, we assume the $x_b$ a priori is accurate enough and only one independent variable $x_m$ is necessary to simplify the process. Since the measurement error of $x_b$ is usually small and unbiased relative to the a priori profiles [23], the results are not sensitive to this simplification. The family of refractivity profiles derived from the Abel-inversion results of Equation (4), with respect to different $x_m$, is shown in Figure 6a. As shown in the figure, the black, dotted curve retrieved by standard Abel-inversion is negatively biased at the ducting altitude compared to the refractivity profile depicted by Equation (1), which is shown as the blue, solid line. The dashed lines in different colors are the family of the refractivity solutions $N_f\left(h,dx\right)$ with different $dx$ values based on the Abel inversion result in Equation (4) and are given ducting height (2 km). As described above, the Abel-inversion process assumes no ducting occurs; therefore, its profile corresponds to the case with $dx=0$. Since all the candidates (including the truth and the Abel-inverse solutions) correspond to the same direct bending angle profile, it is impossible to identify the correct solution in the family solely using direct bending angle observation.

To distinguish each individual profile in the family, the reflected bending angle of each refractivity candidate $N_f$ can be calculated by the forward model of Equation (8). However, two issues have to be addressed in practice. First, RO profiles can be cut-off before the tangent point reaches the surface due to low signal-to-noise ratio (SNR). The RO refractivity profiles do not start from the surface in most of the cases, consequently making the first integral in Equation (8) incomplete. Here, we extrapolate each $N_f$ by fitting the lowest 0.5 km exponentially to extend it to the surface before calculating its reflected bending angle. This is reasonable because it is frequently observed that the moisture within the boundary layer is well-mixed. A simple error analysis regarding the extrapolation is given in Appendix B. Second, the computational cost of the reflected bending forward model at all levels is relatively high. Since the optimization problem of choosing the best candidate profile should be solved iteratively (which will be stated in the next subsection), a more efficient forward model implementation is required. In this research, the approximate approach proposed by Healy and Thepaut [29], which assumes the refractivity profile $N\left(x\right)$ varies exponentially at each interval, is used:

$$\Delta {\alpha }_{R,i}(a)={10}^{-6}\sqrt{2\pi a{p}_i}{N}_{i}exp\{p_i(x_i-a)\}\left[erf\{\sqrt{p_i(x_{i+1}-a)}\}-erf\{\sqrt{p_i(x_i-a)}\}\right],$$

(9)

where i is the level index number, $erf$ is the error function, and $p_i$ is defined as:

$$p_i=\frac{ln(N_i/N_{i+1})}{x_{i+1}-x_i}.$$

(10)

Note that, in the reflected bending angle calculation, the impact parameter a is always below the surface, and unlike the case of the direct bending angle, all levels above the surface are summed. At the height of ducting layer, the variable $p_i$ can fluctuate significantly due to a large gradient of refractivity and small change of the x, which cause the assumption of exponential refractivity curve to be invalid. In this case, when $p_i<0$ (increasing N or decreasing x w.r.t. height h, the latter could occur when $N\left(h\right)$ has a sharp decreasing slope during ducting) or $p_i>2$ (only slight change of x w.r.t. height h at the local maximum of $x\left(h\right)$ when ducting occurs), the linear shape of $N\left(x\right)$ is assumed instead, and Equation (9) should be modified as:

$$\Delta {\alpha }_{R,i}\left(a\right)=-2\sqrt{2a}{10}^{-6}\frac{N_{i+1}-N_i}{x_{i+1}-x_i}(\sqrt{x_{i+1}-a}-\sqrt{x_i-a}),$$

(11)

and the total reflected bending with given impact parameter a below the surface can be calculated as:

$${\alpha }_R\left(a\right)=\sum _i\Delta {\alpha }_{R,i}\left(a\right).$$

(12)

The resulting reflected bending angle of each refractivity profile in the family of Figure 6a is shown in Figure 6b. It is shown in the figure that the reflected bending angle profiles are not the same for each corresponding refractivity profile in the family. Actually, larger refractivity profiles in the family have smaller reflected bending angle profiles at the same level or a larger impact parameter, and vice versa. This characteristic makes the reflected RO signal a suitable constraint to distinguish different profiles below the ducting layer.

Based on Figure 6, the true profile lies between $dx=0.05$ to $dx=0.10$. To search for the exact $dx$ value and its corresponding profile solution, the least square approach is used, as stated in the following subsection.

3.2. Least Squares Estimation

Unlike the approach proposed by Wang et al. [23], which used optimal estimation to estimate both $dx$ and $x_b$ simultaneously, here, $dx$ is the only variable to be estimated which does not have an a priori; hence, a simple non-linear least squares can fit the needs here. The reflected bending angle of the candidate profile ${\alpha }_R^C$ for each iteration is calculated by the forward model described in Section 3.1. The trust region reflective algorithm is used in this research for non-linear least squares implementation to search for the $dx$ that could give the minimum reflected bending gradient difference between the candidate and the observation:

$$dx=argmin{\displaystyle \sum _a}{\left[\frac{d{\alpha }_R^C\left(a\right)}{da}-\frac{d{\alpha }_R^O\left(a\right)}{da}\right]}^2,$$

(13)

where $argmin$ means the argument of the minimum. ${\alpha }_R^O$ is the reflected bending angle observation calculated by PM or other bending angle retrieval methods from the excess phase measurements. The reason that the bending gradient is chosen over the bending angle as the cost function is that, in reality, the refractivity profile cannot be well-represented by the simple bilinear model, and the calculated reflected bending angle, which mostly contributed by ducting when tangent height was located right below the surface, is sensitive to this difference. On the other hand, the bending gradient can be written as [30]:

$$\frac{d{\alpha }_R\left(a\right)}{da}\approx -\sqrt{\frac{1}{2r_E}}\frac{d_{th}}{h_m^{3/2}}+\sqrt{\frac{2}{r_E(a_S-a)}},$$

(14)

where $d_{th}$ is the thickness of the ducting layer. $a_S$ is the impact parameter at the surface, and $r_E$ is the radius of the Earth. The first term described the bending gradient due to atmospheric refraction and is dominated by the ducting structures. But the main contribution of the bending gradient comes from the reflection described by the second term, when tangent heights are located close to (and below) the surface. Therefore, the bending gradient will be insensitive to the difference in ducting structure between model and observation, as well as better-related to the refractivity at surface that is determined by $x_m$. An example using two different criteria will be provided in the next section.

By defining the cost function as Equation (13), the profile with the corresponding reflected bending angle that best matches the slope of observation will be selected. In practice, the reflected RO signal will be significantly attenuated when its frequency deviates from the direct signal because the PRN code no longer aligns with the one locally generated in the open-loop receiver. In a setting case, this happens in the beginning stage of the occultation, which causes the absence of the reflected bending when $a<a_s-0.5$. Additionally, the interference by the direct signal around the surface could spread the spectrum at the impact parameter $a_s$. Therefore, we only compute the bending difference between $a_s-0.4$ km to $a_s-0.1$ km and use it as the observation for least squares.

Figure 7 shows the reconstructed result of the same case stated in the previous section using the least squares approach. The reflected bending angle profiles ${\alpha }_R^O$ (True) and the best candidate ${\alpha }_R^C$ selected by least squares (Fwd) are shown in Figure 7a. The least square chose the $dx$ that leads to the minimum difference of the two curves between the range of $a_s-0.4$ and $a_s-0.1$, which are indicated as black, dashed lines. The resulting refractivity profile is shown in Figure 7b,c. As indicated in the figures, the profiles chosen by least squares are able to reconstruct the multi-valued function $x\left(h\right)$, which was erroneously taken as a single-valued function in Abel inversion. It is shown that the least squares can successfully determine the correct profile in the family and reduce the N-bias due to ducting.

4. Results

4.1. End-to-End Simulation

To verify the reconstruction algorithm described in Section 3, we conducted the same end-to-end experiment as stated in Wang et al. [23] and Xie et al. [18]. The refractivity data set from the radiosonde measurements of Variability of the American Monsoon Systems (VAMOS) Ocean–Cloud–Atmosphere–Land Study (VOCALS) campaign [31], which is located at one of the regions that has the highest frequencies of ducting conditions [32] and leads to a large N-bias in RO refractivity retrievals [33], was used to simulate both direct and reflected bending angles. As shown by Wang et al. [23], 20 cases with strong refractivity gradient at the PBL top have been selected for statistical analysis, and six among them are shown in Figure 8. The blue lines in Figure 8 are refractivity profiles from the radiosonde observations (RAOB) data, which are taken as the truth. It can be observed in these cases that the refractivity profiles has a sharp gradient at the height about 1.5 km. By applying the forward Abel to generate the direct bending angle and simulate the standard RO refractivity observation (black, dashed line) by Abel inversion, the divergence beneath the top of PBL to the surface can be clearly observed. The corresponding direct bending angle of the true profile (RAOB) and the reflected bending angle of both true and Abel inversion profiles calculated by forward Abel are shown in Figure 9.

The reconstructed refractivity for each case is then estimated by the least square approach stated above, based on the negatively biased Abel inverse profiles. This approach looks for the profile candidate with the proxy of $x_m$, in which the reflected bending angle (shown as orange, solid lines in Figure 9) derivative in the shaded area best matches the simulated true reflected bending angle from the RAOB observations (blue, dotted line in Figure 9). The reconstructed refractivity results are shown as orange, solid lines in Figure 8. As the results suggested, the reconstruction can successfully correct the N-bias below the ducting layer compared to the Abel-inverse profile by 5%. Note that the reconstructed profile inside the ducting layer comes from modeling, instead of the data itself, due to the information gap mentioned in the Introduction. Because the assumption of bilinear structure and the constraint of calculating $x_m$ and $h_b$ stated in Appendix A cannot absolutely represent the true thermodynamic structure inside the ducting layer, the resulting profile around the top of PBL (1.5 km) does not align with the RAOB observations. Since the ducting layer only contributes a fraction of the reflected bending angle derivative compared to the surface reflection [30], this discrepancy does not much affect the estimation based on the bending angle derivative constraint, and the reconstructed profile below ducting can still remove most N-bias. However, this discrepancy could affect the result if the bending angle difference, rather than bending angle derivative difference, is used as the cost function. The reconstruction algorithm could erroneously choose the “best-matched” profile, but with a different surface impact parameter, and result in poor estimation. An example of using bending angle difference is provided in Figure 10. It can be observed that the chosen bending angle profile intersects with the simulated observation, which indicates that the gradient is relatively different. The estimation result is thus severely biased down to the surface compared to the true profile.

The statistical reconstruction results compared to the RAOB observations are given in Figure 11. The refractivity difference in percentage in the horizontal axis is defined as:

$$\Delta N\left(h\right)=\frac{N\left(h\right)-N_{RAOB}\left(h\right)}{N_{RAOB}\left(h\right)}\times 100\%,$$

(15)

where $N_{RAOB}\left(h\right)$ is the refractivity observation from RAOB as a function of height h, and $N\left(h\right)$ is either the results of Abel-inversion (black, dotted lines) or the reconstruction (orange, solid lines). As Figure 11 shows, the least squares reconstruction can decrease the N-bias from 5% in Abel inverse profiles to less than 1% at 0.5 km below the top of the ducting layer. The large difference between the reconstruction results and the RAOB observations between $h_t-0.5$ km to $h_t$ is typically due to mis-modeling inside the ducting layer, but this error does not have much influence on the reconstruction below.

4.2. Actual RO Cases

In practice, the reflected bending angle retrieval is not straightforward. The reflected bending angle ${\alpha }_R\left(a\right)$ extracted by PM can be a multi-valued function in the impact parameter domain due to a horizontally non-uniform atmosphere [30]. In addition, the smoothing of the reflected bending angle over the impact parameter in PM can ”shift” the profile vertically because the reflected bending angle has strong gradient close to the surface. To overcome the difficulty, the conventional GO processing for bending angle calculation in the time domain after removing the direct ray path component in the impact parameter domain is a more favorable approach [30]. According to Gorbunov et al. [30], the reflected RO ray is monotonic in the time domain because of the large rate of the grazing angle increasing and is immune to the multi-path effect, which is commonly observed in the direct RO signals over the lower troposphere. Therefore, PM is not essential to the reflected bending angle retrieval in practice, and it is only used in this research as a tool to remove the direct signal in the impact parameter domain. The details of the reflected signal extraction using PM and the application of inverse PM to transform the signal back to time domain for conventional GO processing is described in Appendix C.

The GO algorithm calculating the reflected bending angle is the same as the one used for direct signal described in Gorbunov et al. [30]. The processing is demonstrated by an example with a strong reflection component, as shown in Figure 12a. Compared to its corresponding collocated European Centre for Medium-Range Weather Forecast (ECMWF) model, the RO retrieval shows a large negative bias, which indicates this occultation event is impacted by ducting. By applying the filtering process described in Appendix C, the direct signal is removed and the component left in the spectrogram (Figure 12b) is the extracted reflected signal in the time domain. The resulting amplitude and the excess Doppler example of the reflected signal in the time domain are shown in Figure 12c,d. The blue lines are for full signal, and the orange lines are the reflected signal only. The excess Doppler of full and reflected signals align to each other at the end of setting occultation because the ray path difference is decreasing through the time. Note that, while the reflected SNR reaches 40 V/V in average for this specific case, in most other ducting cases, it is below 30 V/V.

Here, we hand-picked four actual RO cases with a clear reflected signal component in their spectrogram from the Radio Occultation Meteorology Satellite Application Facility (ROM SAF) reflection flag database [34]. All four cases are located in the north-east Pacific Ocean between September and November, where the strong inversion caused stratocumulus deck is frequently observed. The excess Doppler is then used in GO to calculate the reflected bending angle, which is shown in Figure 13 as the blue, dashed line. To suppress the noise, a 1-s window moving average on both the reflected bending angle and its corresponding impact parameter is applied. Only the portion where $0.005rad<{\alpha }_R<0.015rad$ is used as the measurement in this research to avoid the erroneous estimation due to large fluctuation. The best-matched profiles can be found using least squares, shown as the orange lines in Figure 14, and their corresponding reflected bending angle is shown as the orange line in Figure 13. The collocated ECMWF profiles are provided for comparison. It can be observed that the chosen critical height $h_b$ in the reconstructed profiles may be different from the ones that shown in ECMWF profiles, which could be due to the low bias of ECMWF model. But, as shown in the results, the reconstructed profiles significantly reduce the bias observed in original RO profiles (black, dashed line) and closer to the referenced ECMWF (blue, solid line) profiles.

While the unbiased candidate profiles can be identified in these cases, it is not common for the reflected signals to reach such a high SNR. In fact, the resulting reflected bending angle, in most of the cases calculated by GO, contains significant fluctuation, which makes applying least-squares estimation difficult. This large variance of the bending angle is due to the low SNR nature of the reflected signal caused by the delay offset and frequency divergence from the local generated signal, as well as the interference from the direct signal at the end of the setting occultations (or the beginning of the rising occultations). To quantify the uncertainties in the reflected bending calculation, the spectral width in the impact parameter is shown in Figure 13 as the shaded area. The calculation of the uncertainty for this specific case is provided by Gorbunov et al. [35]. It can be observed that the bending angle could contain the impact parameter uncertainty of 0.2 km and noticeably higher when its corresponding impact parameter is close to the surface. This big uncertainty can, therefore, jeopardize the proposed reconstruction method, and possible improvements are proposed in the discussion section.

5. Discussion

As shown in the previous section, using the constraint of the reflected bending angle can effectively reduce the N-bias caused by ducting. While this approach significantly enhanced the usability of the N-bias correction algorithm in practice, some limitations do exist. First, this reconstruction method does not provide additional information of the missing structure that ordinary GNSS-RO cannot detect, such as the profile inside the ducting layer and below the penetration depth. The $x-h$ structure inside the ducting layer and the refrativity profile below the RO cut-off height is provided by the bilinear shape and exponential extension assumption, respectively, rather than retrieved from the RO measurements. Therefore, the error in the reconstruction results can show up due to the erroneous modeling of the thermodynamic structure. While the relationship between them is not straightforward and the sensitivity is unknown, the reflected bending angle can potentially provide the surface refractivity information and the physical constraint inside the ducting layer. Thus, through the better reflected bending angle processing, more information could possibly be extracted and enhance the reconstruction in the future.

Another limitation is the low signal strength of the reflected signal. Although the amplitude information is not directly used in the reconstruction algorithm, low SNR could cause a large variance in the retrieved reflected bending angle. As mentioned above, the SNR of the reflected signals is lower because, in the current open-loop tracking receiver, the local carrier frequency and delay correlator are set to track the direct signal. Therefore, the reflected signal strength decreases when its path length has a larger difference than the direct path length (i.e., higher tangent point altitude). The other possible cause is the signal strength loss due to the sea surface roughness (g), which can be quantified by the Rayleigh roughness parameter [36]:

$$g=4\pi \left(\frac{S_h}{\lambda }\right)sin\varphi ,$$

(16)

where $S_h$ is the standard deviation of the surface height about the local mean, $\lambda$ is the signal wave length, and $\varphi$ is the grazing angle. Equation (16) shows that the small grazing angle significantly decreases g to less than 0.5, which eases the roughness problem. The corresponding reduction factor of the specular component can be modeled [36] as:

$${\rho }_s=exp\left(-1/2g^2\right),$$

(17)

where ${\rho }_s$ is the reduction factor of the reflected signal amplitude. The average global sea surface wind speed is around 6 m/s, which could statistically lead to 1 m significant wave height ($S_h=0.25$) over a fully developed sea [37]. In the reflection case of −300 m tangent point height with grazing angle of 1e-3 $rad$, ${\rho }_s=0.987$ and causes only less than 2% of signal amplitude attenuation. Even on a rough sea surface of 3m significant wave height, ${\rho }_s=0.89$ and still causes only $\sim 1/10$ of amplitude loss. Therefore, only in the most extreme cases could the surface roughness negatively impact the proposed N-bias correction algorithm.

Due to the SNR limitation, the fraction of the high reflected SNR cases among the reflection-detected occultation suitable for the reconstruction algorithm may be considerably less than 45% in mid-latitude regions. This requires further investigations. But the low SNR issue can be overcome in two ways. First, the reflected SNR can be enhanced if the total SNR increases. The SNR of the recently launched COSMIC-2 mission can reach 2500 V/V with the newly designed beam-forming antenna compared to an average of 700 V/V in the current COSMIC mission. Hence, the reflected SNR can be significantly improved, and a higher percentage of RO cases that the proposed method can be applied is expected. The other possible solution is to modify the receiver operation to track the reflected signal with a separate RF channel. The advantage of this approach is that no further processing to extract the reflected signal (Appendix C) is required, which can simplify and improve the reflected signal processing on the ground. These improvements can effectively enhance the reflected SNR and ease the application of the proposed reconstruction method in practice.

6. Conclusions

GNSS-RO is a valuable measurement technique for probing the vertical thermodynamic structure of the PBL from space [6]. However, it is known that a negative refractivity bias (N-bias) can exist within PBL, which will lead to retrieved humidity and temperature biases. Over the stratocumulus regon in the subtropical oceans, the N-bias is dominated by the prevalence of elevated ducts caused by the strong inversion capping the PBL [14,15]. To improve GNSS-RO’s capability of characterizing PBL structure, correcting the N-bias within the PBL caused by ducting is critical. Due to the under-determined condition in standard Abel inversion, the N-bias generated by strong refractivity gradient ($<-157Nunit/km$) can reach 5% within the PBL. In previous studies [18,23], the analytical solutions based on the biased refractivity profile and the physical constraints, including the RO measurements at the surface or external information (e.g., precipitable water) from other collocated observations, are derived. However, actual RO measurements rarely reach the surface, and collocation between GNSS-RO and other observations does not always exist. In addition, the use of external information makes the solution dependent on these data, which may contain biases themselves.

In this article, a new bias reduction method is introduced using the reflected signal component during the occultation events. The grazing components were extracted from the received RO signals, and the bending angle of their ray paths was used as the constraint to solve the ill-posed inversion problem. We implemented the least square algorithm to select the candidate profile that best matches the observed reflected bending gradient from the solution family. The significant advantage of this method is that it does not rely on external information; thus, it eliminates possible errors arising from the assisting observations or models. Based on the simulation and actual RO data results, the proposed reconstruction method is able to reduce the N-bias to less than 1% below the ducting layer. The new construction method has also been applied to the actual RO data: the resulting solutions removed most of the N-bias and are consistent with the collocated model profile. While several limitations, including the low SNR and surface characteristics, could affect the applicability of the proposed method, the challenges can be overcome with the higher SNR and the better grazing RO signal tracking technique in future RO missions.