## Abstract

The global positioning system (GPS) radio occultation (RO) technique has demonstrated the ability to precisely probe earth’s atmosphere globally with high vertical resolution. However, the lowermost troposphere still presents some challenges for the technique. Over moist marine areas, especially in subtropical regions, a very large negative moisture gradient often exists across the thermal inversion capping the marine boundary layer (MBL), which frequently causes superrefraction (SR), or ducting. In the presence of SR, the reconstruction of refractivity from RO data becomes an ill-posed inverse problem. This study shows that one given RO bending angle profile is consistent with a continuum (an infinite number) of refractivity profiles. The standard Abel retrieval gives the minimum refractivity solution of the continuum and thus produces the largest negative bias, consistent with a negative bias often present in the retrieved refractivity profiles in the moist lower troposphere. By applying a simple linear parameterization of the refractivity structure within and just below the SR layer, an analytical relation between the Abel-retrieved refractivity and a continuum of solutions is derived. Combining the Abel retrieval and the analytical relation with some physical constraints, a novel approach is developed to reconstruct the vertical refractivity structure within and below the SR layer. Numerical simulation studies in this paper have demonstrated the great potential of the reconstruction method to provide a much-improved retrieval in the presence of SR, and the method should greatly enhance the ability to measure the MBL structure globally using the GPS RO technique.

## 1. Introduction

The troposphere, where most weather phenomena occur, can be divided vertically into two parts: the planetary boundary layer (PBL) and the free troposphere. The PBL, extending from the surface up to a height range from a few tens of meters to several kilometers, is influenced directly by earth’s surface, responding to forcings such as frictional drag, solar heating, and evapotranspiration. In addition, the PBL controls the evaporation and turbulent redistribution of water into the atmosphere, thus significantly affecting the surface fluxes and the global distributions of cumuliform and stratiform clouds, which in turn influence the radiative forcing of the earth’s climate system (Garratt 1992). Since the PBL is the interval where energy, momentum, and mass are exchanged between the earth’s surface and the free troposphere, it is important to include a realistic representation of the PBL in weather and climate prediction models (Heckley 1985; Albrecht et al. 1986; Betts and Ridgeway 1989).

The PBL over land has been investigated intensively via regular soundings, tower measurements, and various field observations (Lettau and Davidson 1957; Swinbank 1968; Clarke et al. 1971; Izumi 1971; Izumi and Caughey 1976; Clarke and Brook 1979). However, the marine boundary layer (MBL, i.e., PBL over the ocean) is still poorly understood because of the general scarcity of observational data, especially over remote marine regions. The complex physical processes involved in the coupling between the ocean surface and the boundary layer have proven very important and difficult for weather (Hong and Pan 1996; Beljaars and Viterbo 1998) and climate models (Mechoso et al. 1995; Bretherton et al. 2004; Zeng et al. 2004) to simulate.

Global characterization of the MBL requires remote sensing from space. However, due to the relatively short vertical extent of the MBL, and the frequent cloudiness at the MBL top, it is difficult to obtain the MBL structure with most of the current remote sensing techniques. Sensing the MBL via passive infrared (IR) observations is limited due to the inability of IR radiation to penetrate clouds. While nadir-viewing passive microwave sounders can sense the atmosphere over oceans in both clear and cloudy conditions, their vertical resolution is far too poor for boundary layer characterization (Curran 1989). Satellite-borne lidars have used the reflectivity of MBL aerosols to determine the MBL top when the overlying free troposphere is clear (Curran 1989; McCormick et al. 1993; Winker et al. 1996; Cosma-Averseng et al. 2003), but they cannot probe much into the MBL. On the other hand, several features of the global positioning system (GPS) radio occultation (RO) technique suggest that it has a great potential of sensing the MBL (Kursinski et al. 1997; Hajj et al. 2004). These features include global coverage, high vertical resolution, high precision and accuracy, and the ability of GPS signals to penetrate clouds.

While the GPS RO technique is a very promising approach for characterizing the MBL remotely, there are still at least two significant challenges to overcome to enable routine probing of the moist lowermost troposphere. The first is the development of an instrument (receiver) that properly recovers the phase and amplitude of the occulted signal in the presence of large and variable refractivity structure. Doing so requires open-loop GPS receivers with sufficient signal-to-noise ratio (e.g., Sokolovskiy 2001b) analogous to the open-loop receivers used for decades in planetary occultation missions (e.g., Fjeldbo and Eshleman 1968; Eshleman 1973). The second challenge, which is the focus of the present paper, is coupled to the presence of a vertical interval in which the refractivity gradient is so large that the radius of curvature of an internally trapped ray becomes less than the radius of the earth. This is commonly known as superrefraction (SR) or ducting. Under such conditions, there can be no ray tangent point within this altitude interval for an external ray that starts and ends outside the atmosphere (Sokolovskiy 2003), and the standard Abel integral transform does not yield the correct refractivity at this and lower altitudes. In fact, the Abel retrieval will systematically underestimate the actual refractivity profile at and below the height where SR occurs as shown in simulation studies by Sokolovskiy (2003). Superrefraction is expected to occur frequently near the top of the boundary layer over oceans as indicated in numerical weather prediction model analyses (von Engeln and Teixeira 2004) and balloon soundings.

These results are consistent with a negative bias observed in the retrieved refractivity from various GPS RO missions, and noted by a number of researchers beginning with Rocken et al. (1997) who found the bias in the warm and wet regions of the lower troposphere. The negative bias is most severe in the Tropics and at altitudes below 2 km, and can be as large as −10% near the surface (Ao et al. 2003; Beyerle et al. 2003). It is generally believed that the bias arises from a combination of (a) phase-locked loop (PLL) receiver tracking errors that result from the complicated signal dynamics typical for propagation in the moist troposphere (Sokolovskiy 2001a; Ao et al. 2003; Beyerle et al. 2003), and (b) superrefraction (Kursinski et al. 1997, 2000; Sokolovskiy 2003).

While the solution to the PLL receiver tracking errors appears largely in hand with the open-loop GPS tracking being implemented at the Jet Propulsion Laboratory (JPL) for the Constellation Observing System for Meteorology Ionosphere and Climate (COSMIC) mission (Rocken et al. 2000), the superrefraction problem is still in need of a solution. Here we present an approach to the SR problem that offers to solve what has been a fundamental retrieval problem in remotely sensing the moist boundary layer using the GPS RO technique.

There are six parts in this paper. In section 2, the effect of SR in geometrical optics and the resulting negative refractivity bias are addressed. In section 3, a novel reconstruction method is proposed to solve the SR and so-called *N*-bias problem. In section 4, simulations are used to demonstrate the reconstruction concept, and preliminary error analysis results are presented. In section 5 we discuss remaining outstanding problems related to the reconstruction method and suggest directions for future work. We summarize and conclude in section 6.

## 2. Refractivity bias due to superrefraction

### a. Superrefraction in the marine boundary layer

In the neutral atmosphere, the refractivity *N* (measured in dimensionless so-called *N* units) at GPS frequencies [*N* = (*n* − 1) × 10^{6}, where *n* is the refractive index] is related to the atmospheric pressure (*P* in mbars), temperature (*T* in kelvins), and water vapor partial pressure (*P _{w}* in mbars) through (Smith and Weintraub 1953)

Over moist marine areas, especially in subtropical regions, a very sharp temperature inversion and a large negative moisture gradient are often observed at the top of the MBL. These sharp temperature and moisture gradients give rise to a large negative refractivity gradient that causes a large bending of GPS RO signal paths. It is called superrefraction, when the vertical refractivity gradient exceeds a critical value, that is, *dN*/*dz* < −1/*r _{E}* ≈ −157

*N*-unit km

^{−1}(

*r*is the curvature radius of the earth). At the top and the bottom of the SR layer, where the refractivity gradient just equals the critical value (i.e., critical refraction), the ray path in principle becomes circular and propagates at a fixed height above the earth’s surface.

_{E}Figure 1 shows a representative high-resolution dropsonde sounding off the coast of southern California in July 2001, during the Dynamics and Chemistry of Marine Stratocumulus Phase II (DYCOMS-II) boundary layer observational campaign. It shows a well-defined boundary layer capped by a strong temperature inversion (Fig. 1a) and sharp negative moisture gradient (Fig. 1b) between 1 and 1.2 km above the ocean surface. According to (1), the resulting refractivity (Fig. 1c) and its gradient (Fig. 1d) profiles are shown. Note that the refractivity profile is smoothed to 50-m vertical resolution and thus the major refractivity gradient at the MBL top can be easily distinguished. The large vertical moisture and temperature gradients lead to a sharp refractivity gradient to which RO signals are very sensitive. Compared to the temperature changes, the specific humidity profile exhibits a larger fractional change across the MBL top, and typically dominates the refractivity behavior in this interface. Note that the maximal negative refractivity gradient greatly exceeds the critical value, which indicates strong SR at the top of this MBL.

### b. Standard retrieval approach in the geometrical optics approximation

A schematic plot of the GPS RO geometry is shown in Fig. 2. Under the assumption of spherical symmetry, a ray satisfies Bouger’s (or Snell’s) law (Born and Wolf 1964) in the geometrical optics approximation, that is,

where *r* is the distance from the center of curvature, *ϕ* is the angle between the ray path and the radial direction, and the impact parameter *a* is a constant for a given ray. For a given ray with a well-defined impact parameter, the total refractive bending angle *α* as a function of *r _{t}* (i.e., the curvature radius of the atmosphere at the tangent point), is given as (Fjeldbo et al. 1971)

Normally, by using the substitution *x* = *n*(*r*)*r* and *a* = *n*(*r _{t}*)

*r*, (3) can be transformed to (4) when

_{t}*x*(

*r*) is a monotonic function:

The refractive index profile *n*(*x*) can be solved for from (4) via substitution of variables and the Abel transform (Fjeldbo et al. 1971):

In the case of SR, *x*(*r*) is not a monotonic function inside the SR layer and a layer right below it (called the shadow layer). Thus, (4) will be invalid inside and below the SR layer (Sokolovskiy 2003). However, using the true refractivity profile (solid line in Fig. 1c) as an input, we can still apply (4) above the SR layer, whereas for the layer below the shadow layer, we will need to use (3) instead of (4) to calculate *α*(*a*). Having the *α*(*a*), we can apply (5) to obtain the Abel-retrieved refractivity profile (dashed line in Fig. 1c). The Abel-retrieved refractivity profile is identical to the true profile above the higher critical refraction point, which marks the upper limit of the SR layer. However, the SR results in a negative bias in the Abel-retrieved *N*(*r*) within and below the SR layer. Therefore, SR represents a fundamental limitation for the standard Abel inversion.

## 3. Reconstruction method

Sokolovskiy (2003) pointed out that the Abel retrieval and the true refractivity profiles both yield the same bending angle profile in the case of SR. Our further research reveals that when SR occurs, there is a continuum (an infinite number) of refractivity profiles that yield the same bending angle profile. The standard Abel solution is the minimum refractivity solution in the continuum and therefore produces the largest negative bias. Thus, reconstructing the refractivity profile from the bending angle in the presence of SR is an ill-posed inverse problem.

### a. Analytical relation between the Abel retrieval and a continuum of refractivity profiles

As already discussed in section 2, SR is caused primarily by the large negative moisture gradient across the thermal inversion capping the MBL. Analysis of many high-resolution sounding profiles over oceans indicates that, in most cases, there is one and only one major SR layer that exists at the top of the MBL. In Fig. 3, a schematic MBL structure with an elevated SR layer is shown. The SR layer has two critical refraction points at *h*_{2} and *h*_{3}, or *r*_{2} and *r*_{3} (tangent height: *h* = *r* − *r _{E}*). We separate the portion of the true refractive index profile below

*h*

_{3}into three intervals:

A, the mixing layer: surface to

*h*_{1}in which*h*and*n*are referred to as*h*_{A}and*n*_{A},B, the shadow layer:

*h*_{1}to*h*_{2}in which*h*and*n*are referred to as*h*_{B}and*n*_{B},C, the SR layer:

*h*_{2}to*h*_{3}in which*h*and*n*are referred to as*h*_{C}and*n*_{C}.

As seen in Fig. 3a, within the shadow layer and the SR layer (*h*_{1}–*h*_{3}), *h*(x) is a multivalued function. As the tangent points of rays descends toward the height of critical refraction (i.e., *h*_{3}), the bending angle becomes infinite (Fig. 3c). Moreover, there are no external rays that can have tangent points inside the SR layer or the shadow layer (Sokolovskiy 2003). Note that the impact parameter at *h*_{1} and *h*_{3} are both equal to *x*_{1}, and the SR and shadow layer collapse to a singularity (i.e., an infinite large bending angle) at *x*_{1} in the bending angle profile (Kursinski et al. 2000). The Abel-retrieved refractive index profile is indicated by dashed lines, that is, *h̃*(*x*), *ñ*(*r*), and *ñ*(*x*) in Figs. 3a, 3b and 3d, respectively.

According to Sokolovskiy (2003), both the true and the Abel-retrieved refractivity profiles yield the same bending angle profile. Therefore, according to (3), for a given impact parameter, *a* < *x*_{1}, we have

where *a* = *ñ*(*r̃ _{t}*)

*r̃*=

_{t}*n*

_{A}(

*r*)

_{t}*r*. The left-hand-side (lhs) is for the Abel-retrieved profile and the right-hand-side (rhs) is for the true profile. The lower limit

_{t}*r̃*

_{t}of the Abel retrieval refers to a height below

*r*

_{3}, whereas

*r*of the true profile refers to a height below

_{t}*r*

_{1}. The limits,

*r̃*and

_{t}*r*, must correspond to the same impact parameter to produce the same

_{t}*α*(

*a*) profile. Since the true and the Abel-retrieved refractive index profiles in Fig. 3b are identical above

*r*

_{3}(or equivalently

*h*

_{3}), the corresponding parts of the integrals in (6) above

*r*

_{3}are identical. Thus, we need only consider the parts of the integrals below

*r*

_{3}, and instead of (6), we can write

Here we have broken up the rhs integral into three parts, one for each of the intervals A, B, and C (cf. Fig. 3). Rewriting (7) in terms of *x* (using the substitution *x* = *nr*), we get

Next, by taking advantage of the fact that *x* and *h* uniquely define *n* [i.e., *x* = *nr* = *n*(*h* + *r _{E}*)], we can simplify (8) in terms of height

*h*instead of refractive index. We also make an important assumption/parameterization that the dependence of height on

*x*across the SR and shadow layers [i.e.,

*h*

_{C}(

*x*) and

*h*

_{B}(

*x*)] can be approximated by two straight-line segments:

This simplification is based on the examination of many high-resolution radiosonde and dropsonde soundings over oceans, and yields a relatively simple analytical solution for (8). In terms of *h* versus *x*, the solution can be written as (see the appendix for details)

where *z* = *x*_{1} − *x*/*x*_{2} − *x*_{1}.

Given *x*_{1}, *x*_{2}, *h*_{1}, and *h*_{2}, (10) relates the Abel-retrieved refractivity profile *h̃*(*x*) to a profile below *h*_{1}, that is, *h*_{A}(*x*). This implies that the Abel-retrieved profile does not contain explicit information about the vertical structure inside the SR and shadow layers [i.e., *h*_{B}(*x*) and *h*_{C}(*x*)]. This makes sense, since the Abel-retrieved profile is derived from the bending angle profile, which only contains implicit information about the vertical structure inside the SR and shadow layers (note that rays with tangent point below *h*_{1} travel through these layers, and the bending angle below *x*_{1} therefore depends on the vertical structure inside the SR and shadow layers to some extent).

In the presence of SR, a large positive spike in the bending angle and the corresponding large negative gradient (close to critical) in the Abel-retrieved refractivity can be used as a potential indicator of SR (Sokolovskiy 2003). On the other hand, an abrupt weakening (or disappearance) of the RO signal intensity at the receiver potentially provides the information of the existence of the SR layer (Hajj et al. 2004). In principle, the impact parameter (*x*_{1}) and the upper limit of the SR layer (*h*_{3}) can be obtained from the bending angle profile and the corresponding Abel-retrieved refractivity profile, respectively. However, for real occultation data, when the bending angle contains errors, identifying *h*_{3} and *x*_{1} could become a challenging problem and will be discussed in section 5. For simplicity, in the following we assume *h*_{3} and *x*_{1} are known parameters obtained from the bending angle and refractivity retrievals.

Knowing *h*_{3} and *x*_{1}, (10) has two degrees of freedom, *h*_{1} and *x*_{2}. Based on the simple linear segment assumption for *h*_{B}(*x*) and *h*_{C}(*x*), there is one more degree of freedom left for solving *h*_{B}(*x*) and *h*_{C}(*x*), namely, the lower limit of the SR layer, *h*_{2}. To emphasize the degrees of freedom we write (10), (9a), and (9b) as follows:

Besides the independent variable, *x*, *h*_{B}(*x*), and *h*_{C}(*x*) depend on *x*_{2}, *h*_{1}, and also *h*_{2}, whereas *h*_{A}(*x*) does not depend on *h*_{2} but only depends on *x*_{2} and *h*_{1}. Equations (11a)–(11c) relate the Abel-retrieved profile to a continuum (an infinite number) of solutions. All these solutions correspond to the same bending angle profile and can be described by *x*_{2}, *h*_{1}, and *h*_{2}. However, not all combinations of *x*_{2}, *h*_{1}, and *h*_{2} make sense physically. For example, it is obvious from the formulation of the problem that *x*_{2} ≥ *x*_{1}, and that *h*_{1} ≤ *h*_{2} ≤ *h*_{3}. It is interesting to point out that the second term on the rhs of (10), or equivalently *G*_{h1, x2}(*x*) in (11a), actually is a measure of the error in the Abel-retrieved refractivity profile, when *h*_{A}(*x*) represents the “truth.” It is easily verified that *G*_{h1, x2}(*x*) (*x*_{0} ≤ *x* ≤ *x*_{1}) is always nonpositive, and so *h̃*(*x*) ≥ *h*_{A}(*x*) and *ñ*(*x*) ≤ *n*_{A}(*x*). Therefore, the Abel solution is the minimum refractivity profile in the continuum of solutions and thus produces the maximum negative refractivity error. The error of the Abel solution is affected by three factors, such as (*h*_{3} − *h*_{1}), (*x*_{2} − *x*_{1}), and (*x*_{1} − *x*), which more or less represents the total thickness of the SR and shadow layers, the “strength” of the refractivity gradient inside the SR layer, and the distance below the height of the critical refraction, respectively. Clearly, the error in the Abel-retrieved refractivity profile is proportional to the thickness of the SR and shadow layer (*h*_{3} − *h*_{1}). Also note that the function −*G*_{h1, x2}(*x*) increases monotonically as *z* decreases. Accordingly, the larger the refractivity gradient inside the SR (i.e., *x*_{2} − *x*_{1}), the larger the error is, but generally the error decreases at lower altitudes (i.e., smaller *x*).

With the three degrees of freedom in (11a)–(11c), we need additional constraints to estimate the refractivity structure inside and below the SR layer. There are two steps left for our reconstruction method: (a) derive a subset of a continuum of solutions that are consistent with a given bending angle profile and (b) isolate the solution closest to the “truth.”

### b. Identifying the solution closest to the “truth” from a continuum of solutions

According to many high vertical-resolution sounding observations over global oceans, we find that it is reasonable to make the following assumption, which will act as a constraint in the reconstruction method: the slope of *h*_{A}(*x*) is finite and continuous at the bottom of the shadow layer, that is,

In addition, we shall assume that the bending angle profile retrieved from the RO data extends all the way to the earth’s surface, such that the minimum impact parameter *x*_{0} from the retrieved bending angle profile can be used as a surface constraint. Thus, as a second constraint we have

For real occultation data, however, tracking might stop before reaching the surface, or the impact parameter in the retrieved bending angle profile corresponding to a tangent point near the earth’s surface might not be sufficiently well determined to act as a substantial constraint. Additionally, surface reflections and edge diffraction might also result in an erroneous estimate of *h̃*(*x*_{0}). Possible alternatives to the constraint in (13) will be discussed in section 5.

Equation (12) indicates that *h*_{B}(*x*) ≈ *h*_{A}(*x*) = *h*_{1} − (*h*_{2} − *h*_{1})*z*^{2} for small *z* (*x* close to *x*_{1}). Therefore, Taylor expansion of (10) yields the following equation valid for small *z*:

Invoking the definition of *z* = *x*_{1} − *x*/*x*_{2} − *x*_{1}, we see that to the first order, *h̃*(*x*) − *h*_{3} behaves like a square root function in *x* near *x*_{1}. It follows that

The constant *C*_{0} can be computed from the lhs of (15). Thus (15) potentially gives another relation between *x*_{2} and *h*_{1}, which can be used as a third constraint for the reconstruction method. In general, for a given *C*_{0}, (13) and (15) constitute two nonlinear equations that can be solved simultaneously for *x*_{2} and *h*_{1}. Subsequently, *h*_{A}(*x*) can be found via (10). However, *C*_{0} has to be known quite accurately if we want the square root term in *h̃*(*x*) near *x*_{1} to cancel with a similar term in *G*_{h1, x2}(*x*) such that *h*_{A}(*x*) ends up having finite derivative at *x* = *x*_{1}. In fact, the accuracy of *C*_{0} calculated from the lhs of (15) appears to be insufficient. Numerically, we found that a small positive error in *C*_{0} results in a sharp decrease of the slope of *h*_{A}(*x*) near *x*_{1}, while a small negative error in *C*_{0} leads to a sharp increase of the slope. Consequently, instead of using the lhs of (15), we determine *C*_{0} as the value for which linear regression of the corresponding *h*_{A}(*x*) result in the smallest root-mean-square (RMS) residual. In the present work, the linear regression of *h*_{A}(*x*) was performed over 200 m below *x*_{1}. After *C*_{0} is estimated, the constraint (13) can be relaxed, and a whole family of (*h*_{1}, *x*_{2}) pairs can be obtained from (15) under the constraint *h*_{A}(*x*_{0}) ≤ *h̃*(*x*_{0}). It is worth noting that *C*_{0} is estimated only once using the constraint *h*_{A}(*x*_{0}) = 0.

Each set of *h*_{1} and *x*_{2} thus solved for uniquely defines the structure of *h*_{A}(*x*) according to (10). Because of the straight-line parameterization of *h*_{B}(*x*) and *h*_{C}(*x*), and the continuous slope constraint (11), we can determine *h*_{2} simply by linear extrapolation of *h*_{A}(*x*) until *x* = *x*_{2}. Thus *h*_{B}(*x*) and *h*_{C}(*x*) are obtained via (9a) and (9b). For different values of *h*_{A}(*x*_{0}), this process builds up a subset of a continuum of profiles, which yield the same bending angle profile. In the ideal case of spherical symmetry and perfect receiver tracking, the signal ray path having the minimum impact parameter *x*_{0} would be the lowest ray path with the tangent point touching the earth’s surface. Thus, the solution satisfying *h*_{A}(*x*_{0}) = 0 would presumably be the one in the continuum of profiles that best represents the vertical structure within the MBL. Again, however, we stress that the *h*_{A}(*x*_{0}) = 0 constraint is unlikely to be valid when it comes to real occultation data and is only used in this paper to demonstrate the reconstruction approach in its simplest form.

## 4. Simulation results

### a. End-to-end simulation system

Combining the forward integration and the standard Abel inversion with the reconstruction method, we perform an end-to-end simulation experiment based on geometrical optics as shown in Fig. 4. The end-to-end simulation system includes three parts: (a) a forward part that generates the bending angle profile through the forward integration (3) of an input refractivity profile; (b) an inverse part that retrieves a refractivity profile via Abel inversion using (5); and (c) a reconstruction part that derives a subset of a continuum of refractivity profiles corresponding to the bending angle profile, and isolates one solution that is presumably closest to the truth (the input profile). By comparing the retrieved refractivity profiles with the input profile, we quantitatively evaluate the refractivity errors due to the SR in the Abel inversion and the reconstruction method.

For illustration of the end-to-end simulation system, we use the dropsonde refractivity profile shown in Fig. 1c as an input. The height of the upper limit of the SR layer *h*_{3}, and its corresponding impact parameter *x*_{1}, is obtained from the input *x*(*r*) profile. Note that we smooth the input refractivity profile by applying a 50-m smoothing window. The smoothing removes one thin SR layer (with only about 20-m thickness) above the major SR layer. As seen in many high-resolution radiosonde profiles, some local effects such as turbulence or advection could be the reasons for such thin SR layers, and we believe that they are unlikely to extend horizontally over a large scale, and might not be detectable by the GPS RO technique.

Figure 5 demonstrates the whole concept of the reconstruction method. After applying the forward integration followed by the Abel inversion, the standard Abel-retrieved refractivity profile is achieved (dashed line in Figs. 5a and 5c). Then a subset of a continuum of solutions (dotted) that are consistent with the same bending angle profile (Fig. 5d) can be obtained by the reconstruction method. Finally, the one solution (dash–dot) that has the lowest tangent point touching the earth’s surface [i.e., *h*_{A}(*x*_{0}) = 0] is the solution that is picked out as the reconstructed profile. The reconstructed profile [in terms of *h*(*x*) or *N*(*h*)] shows a very close match with the input or the “true” profile (solid), with only a small deviation within the SR and shadow layers. Figure 5b shows the refractivity errors for both the Abel and the reconstructed solutions. The results demonstrate that the reconstruction method greatly reduces the negative refractivity errors in the Abel-retrieved profile from the maximum of about −8% to less than −1.5% overall. Moreover, the reconstructed profile shows almost a perfect match with the input profile below the shadow layer (around 0.8 km), indicating that the reconstruction method has the potential to eliminate the *N* bias associated with SR.

### b. Case studies

The dropsonde profile collected off the coast of southern California (Fig. 1) is a typical example of MBL structures with an elevated SR layer. This type of MBL is frequently capped by a sharp temperature inversion layer and often has persistent stratocumulus clouds below the inversion layer. The stratocumulus clouds favor the eastern subtropical oceans, where upwelling and cold currents keep the surface waters cool compared to the warm, dry subsiding air aloft (Klein and Hartmann 1993). A few more typical MBL structures with one SR layer from other places are shown in Figs. 6 and 7. These high-resolution observational data have been selected from various field campaigns and island soundings over global oceans. In Figs. 6 and 7, each row demonstrates a simulation case for one sounding profile. The parameter *x* [=*n*(*r*)*r*] versus height, refractivity versus height, and the percentage errors for both the Abel and reconstructed profiles are shown in the three panels from left to right, respectively.

Figure 6 shows four different types of MBL structures having different thickness of SR and shadow layers, and different strength of refractivity gradients in the SR layer. As pointed out in section 3a, these different MBL structures will result in a rather large variation in the magnitude of the Abel-retrieved refractivity errors. The errors of the Abel retrieval show the maximum negative value ranging from less than −5% to over −10% inside the SR or shadow layers, while decreasing rapidly at lower altitudes. However, the errors are significantly reduced after reconstruction for all different types of MBL. The detailed refractivity structures due to the complicated moisture behavior are recovered quite well, and almost no refractivity bias exists in the reconstructed profiles below the shadow layer.

Overlaying the cool southeast Pacific Ocean is the most persistent subtropical stratocumulus deck in the world (Bretherton et al. 2004). The case shown in Fig. 6a is observed in this region off the coast of Peru in South America (11.88°S, 88.07°W). The MBL in this region is similar to what is observed off the coast of southern California. The persistent SR layer in these regions (von Engeln and Teixeira 2004) normally has a very large negative refractivity gradient across a sharp inversion layer at the MBL top. The simulated Abel-retrieved refractivity profile indicates a large negative bias inside and below the SR layer due to the sharp refractivity gradient at the MBL top. In this specific case, the refractivity has the minimum error (−4%) close to the surface, and reaches the maximum error (about −12%) inside the SR layer. Figure 6b shows an island sounding from the HILO site (19.72°N, 155.07°W) in Hawaii. It is over the relatively warm subtropical ocean surface where subsidence in the free troposphere takes place almost all year round. The SR is also very persistent in this region. However, compared to the MBL structures in stratocumulus regions, Hawaiian soundings show relatively higher SR layers, but have less pronounced inversion across the MBL top. Figure 6c shows a radiosonde sounding from Yap Island (9.48°N, 138.08°E) in the tropical western Pacific Ocean. The abundant water vapor available near the warm ocean surface and the active turbulent mixing or convections make the moisture structure in the MBL very complicated, which can be seen in the refractivity profile. Figure 6d shows a representative MBL sounding from St. Paul Island (57.15°N, 170.22°W) in the North Pacific Ocean (close to the Arctic Ocean). The sea surface temperature is relatively cold, and the moisture structure shows much quieter behavior in this region as compared with the subtropical and tropical marine soundings. As a result, the MBL in this region shows thinner SR and shadow layers with a smaller refractivity gradient, and thus a relatively smaller negative *N* bias in the Abel-retrieved profile.

Note again, two assumptions have been made in our reconstruction method: (a) straight-line segment assumption of *h*(*x*) in both SR and shadow layers; (b) continuous slope of *h*(*x*) at the bottom of the shadow layer. Two cases where these assumptions do not hold well are shown in Fig. 7. Figure 7a (33.08°N, 16.35°W) shows a cloud-topped MBL with nonlinear structure inside the SR layer, which is caused by the complex moisture structures inside a cloud. Figure 7b shows another case (31.25°N, 121.75°W) that does not follow the continuous slope assumption very well. Nevertheless, the errors of the Abel retrieval are generally reduced in both cases, and the bias below the shadow layer is virtually eliminated.

In stratocumulus regions (e.g., Fig. 7b), a cloud-topped MBL with a well-mixed layer below the cloud base is often observed. In such situations, the water vapor remains unsaturated below the cloud layer and becomes saturated inside clouds. Therefore, a slope discontinuity in the specific humidity (or relative humidity), and thus in the refractivity, may be present at the cloud base. The discontinuous slope at the cloud base may become an issue when the cloud base happens to be at the same height or above the bottom of the shadow layer. In these cases, extra positive errors may be present inside the shadow layer of the reconstructed refractivity. Simultaneously, a negative error may show up in the estimation of *h*_{2}, which can result in relatively large negative errors inside the SR layer. Figure 7b is a case where the cloud base happens to be at approximately the same height as the bottom of the shadow layer.

### c. Error analysis

To better illustrate how well the reconstruction method works in general, we have selected 55 high-resolution marine sounding profiles smoothed to 50-m vertical resolution. By running these sounding profiles through the end-to-end simulation system, we obtain error estimates for both the Abel and the reconstruction solutions, as shown in Fig. 8. The height of the MBL (normally coincident with the upper limit of the SR layer) varies from about 500 m to near 4 km and the thickness of the SR and shadow layers varies from a few tens of meters up to about 400 m. To make the errors for different cases comparable, we have offset the error plots by the bottom height of the shadow layer (i.e., *h*_{1}).

As seen from Fig. 8, the errors of the Abel retrieval and the reconstruction errors show very different patterns. The errors in the Abel retrieval are always negative, while the reconstruction errors show both positive and negative values. The maximum negative errors for both Abel and reconstructed profiles occur inside the SR layer. The negative errors in the Abel retrieval decline from the maximum and reach the minimum near the surface. The reconstruction errors, however, often reach a maximum positive error inside the shadow layer and then decrease to very small values below the shadow layer. This reconstruction error pattern is primarily caused by the simple linear parameterization for the SR and shadow layers, which result in a sharp kink at *h*_{2}. However, the dropsonde and radiosonde profiles normally bear a rounding edge at *h*_{2}. Thus, the reconstruction method tends to overestimate the refractivity below *h*_{2} and underestimate it inside the SR layer. Also note that, due to the large refractivity gradient inside the SR layer, a small error in *h*_{2} could result in a relatively large error within the SR layer.

Despite of the relatively simple linear parameterization, the reconstruction method gives significantly smaller errors than the Abel retrieval inside the SR and shadow layers. More importantly, much better reconstruction results are achieved below the shadow layer. The errors shrink from the maximum of about −9% in the Abel-retrieved profiles to about ±1% or less near the bottom of the shadow layer in the reconstructed profiles and decrease further below.

## 5. Discussion

The simulation results in this paper demonstrate that the reconstruction method has the potential to eliminate the negative refractivity bias due to SR in the standard Abel retrieval. However, before the approach can be applied to real occultation data, many outstanding problems need to be solved. This section is dedicated to a discussion of these problems.

The sounding profiles with an elevated SR layer show very different MBL structures over different areas of global oceans. Simulation studies on the selected high-resolution marine sounding profiles indicate a large variation of the negative bias in the Abel-retrieved refractivity profile, which can be as large as about −15% in certain situations (Fig. 8). The spatial and temporal variations of the SR structures over oceans based on the sounding profiles need to be investigated. This would be helpful for the assessment of the possible bias in the Abel retrieval and the verification of the assumptions and constraints used in the reconstruction method. It would also help identify the existence of SR in real occultation data.

The bending angle profiles used in this paper are generated by the ideal forward Abel integration. By replacing the forward integration part in the end-to-end simulation system (Fig. 4) with the diffraction-corrected wave-optics simulation followed by the canonical transform (Gorbunov 2002) or the full spectrum inversion (Jensen et al. 2003), we could generate more realistic bending angle profiles with various factors, such as receiver tracking and large-scale inhomogeneous atmospheric effects considered. Then, analyzing the reconstruction results based on these bending angle profiles could lead to the development of a more robust reconstruction method.

We have assumed that the impact parameter (*x*_{1}) and the upper limit of the SR layer (*h*_{3}) can be obtained from the bending angle and refractivity retrievals, respectively. In real occultation data, however, imperfect receiver tracking, the retrieval algorithms (e.g., canonical transform or full spectrum inversion) and large-scale atmospheric horizontal gradients may lead to difficulties of identifying the existence and the height of the critical refraction point directly from the bending angle or signal intensity data. Thus, the errors in the retrieved bending angle, contributed by different sources, need to be examined and separated. Sensitivity analysis of the reconstruction method due to errors in the estimated parameters and assumptions should be investigated.

A robust method to detect the existence of an SR layer and its upper boundary height must be developed before the reconstruction method can be applied to real occultation data. Note that the open-loop receiver tracking technique (Sokolovskiy 2001b) is anticipated to significantly reduce the receiver tracking errors. Thus, with reduced receiver tracking errors, a relatively large bending angle or a large negative refractivity gradient (close to the critical value) in the lower troposphere, in combination with a large negative refractivity bias compared with a numerical weather prediction (NWP) model analysis close to the surface will be a good indicator of the existence of SR. In addition, climatological information of SR over global oceans from sounding data or model analyses could also help identify the existence of SR. With the confidence of the existence of SR, we may then assume that the height of SR (*x*_{1} or *h*_{3}) is close to the height of the maximum bending angle.

Because of the errors in the retrieved bending angle, the lowest impact parameter from the retrieved bending angle profile might not act as a substantial constraint for the reconstruction method. Instead, it may be possible to use an NWP model analysis such as the National Centers for Environmental Prediction (NCEP) or the European Centre for Medium-Range Weather Forecasts (ECMWF) analyses, or use a simple atmospheric surface layer model to infer the near-surface atmospheric structure, and thus the refractivity and the impact parameter close to the surface. Another potential constraint for the reconstruction method, which would be less sensitive to the surface air conditions, could include the integral of water vapor or refractivity inside the MBL from the NWP model analysis.

## 6. Conclusions

The GPS RO technique has the potential to reveal high-resolution vertical structure in the MBL across and below an SR layer. However, in the presence of SR, the problem is ill-posed, and the standard Abel retrieval results in a nonunique solution. That is, one given bending angle profile is consistent with a continuum (an infinite number) of refractivity profiles. The standard Abel transform yields the minimum refractivity solution in the continuum and thus produces the largest negative bias inside and below the SR layer.

Based on a simple linear parameterization of the refractivity structure within the SR and shadow layers, we have developed an analytical relation between the Abel-retrieved refractivity and a continuum of solutions, which are all consistent with one bending angle profile. In combination with a few additional constraints, a subset of the continuum of refractivity profiles can be reconstructed and one solution closest to the “truth” can be isolated.

The error analysis in this paper indicates that the reconstruction method has the potential to provide unbiased refractivity solutions from RO data in the presence of SR, at least below the shadow layer. Negative refractivity errors in the Abel-retrieved profiles are significantly reduced after applying the reconstruction method in an end-to-end simulation experiment. The remaining errors inside the SR and shadow layers are mainly due to the simple linear parameterization used in the reconstruction method.

We have shown very encouraging reconstruction results in this paper, which are based on the examination of bending angle data generated from the forward integration of an input refractivity profile. However, before the reconstruction method can be applied to real occultation data, some of the outstanding issues discussed in section 5 should be addressed.

In conclusion, the proposed reconstruction method should greatly enhance our ability to measure the MBL structure globally using GPS RO techniques. Furthermore, it could lead to an improvement upon the MBL parameterizations used in weather and climate models and eventually help the verification and prediction for global models.

## Acknowledgments

This research was supported by NPOESS and later by NSF. Dr. Bjorn Stevens and Ms. Verica Savic-Jovcic at UCLA are thanked for providing DYCOMS-II soundings. Mr. Michael Brunke is thanked for providing the EPIC dataset. Some ASTEX sounding profiles were obtained from the NASA Langley Research Center Atmospheric Sciences Data Center. TEPPS data were received from Mesoscale Group at the University of Washington. Some island soundings were acquired from the SPARC data center at the Stony Brook University (SUNY).

## REFERENCES

**,**

**.**

**,**

**.**

**,**

**.**

**,**

**,**

**,**

**,**

**.**

**.**

**,**

**,**

**.**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**.**

**,**

**.**

**,**

**,**

### APPENDIX

#### Analytical Relation between the Abel Retrieval and a Continuum of Refractivity Profiles

We start by rearranging and simplifying (8) as follows:

Since *n*(*x*) = *x*/[*r _{E}* +

*h*(

*x*)] and

*ñ*(

*x*) =

*x*/[

*r*+

_{E}*h̃*(

*x*)], and because

*h*≪

*r*, (A1) can be approximated as

_{E}Using *y* instead of *a*, and writing Δ*h*(*x*) = *h*_{C}(*x*) − *h*_{B}(*x*), we define the function *f* (*y*) as

We then perform the following integral transformation (similar to the Abel transform):

Having *h*_{A}(*x*_{1}) = *h*_{1} and *h̃*(*x*_{1}) = *h*_{3} (cf. Fig. 3a), we arrive at

Equation (A4) is a general expression relating *h*_{A}(*x*′) to *h̃*(*x*′), given *h*_{1}, *h*_{3}, and the function *f* (*y*), which depend on the vertical structure inside the SR layer and the shadow layer.

From the definition of *f* (*y*) in (A3) we have

Integrating the last expression by parts yields

In the following we find an approximate solution to (A5) using the assumption that Δ*h*(*x*) can be approximated by two straight lines as described by (9a) and (9b) in the main text. Thus, we have Δ*h* = *h*_{C} − *h*_{B} = *A*(*x*_{2} − *x*), where *A* = (*h*_{3} − *h*_{1})(*x*_{2} − *x*_{1})^{−1}. We may now write (A5) as

With *t* = *x*^{2} − *x*^{2}_{1}, the second term on the rhs of (A6) can be written as

Using the substitution *t* = tan^{−1}(*x*^{2} − *x*^{2}_{1}/*x*^{2}_{1} − *x*′^{2}), the third term on the rhs of (A6) can be written as

Since *x*^{2}_{1} − *x*′^{2} ≪ *x*^{2}_{1}, we Taylor expand the square root in the integrand on the rhs of (A8) to first order and find a relatively simple analytical solution. We get

To be consistent with the approximation made in getting from (A8) to (A9), this is further reduced to

where *z* = *x*_{1} − *x*′/*x*_{2} − *x*_{1}.

## Footnotes

* Current affiliation: Department of Earth and Atmospheric Sciences, Purdue University, West Lafayette, Indiana

*Corresponding author address:* Feiqin Xie, Department of Earth and Atmospheric Sciences, Purdue University, 550 Stadium Mall Drive, West Lafayette, IN 47907-2051. Email: xief@purdue.edu