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) × 106, where n is the refractive index] is related to the atmospheric pressure (P in mbars), temperature (T in kelvins), and water vapor partial pressure (Pw in mbars) through (Smith and Weintraub 1953)

 
formula

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/rE ≈ −157 N-unit km−1 (rE 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.

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.

Fig. 1.

(a)–(d) Dropsonde sounding acquired off the coast of southern California (31.25°N, 121.75°W) at 1350 ∼ 1404 UTC 17 Jul 2001. Thin vertical dashed line in (d) indicates critical refraction.

Fig. 1.

(a)–(d) Dropsonde sounding acquired off the coast of southern California (31.25°N, 121.75°W) at 1350 ∼ 1404 UTC 17 Jul 2001. Thin vertical dashed line in (d) indicates critical refraction.

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,

 
formula

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 rt (i.e., the curvature radius of the atmosphere at the tangent point), is given as (Fjeldbo et al. 1971)

 
formula
Fig. 2.

Illustration of the GPS RO geometry. The radius of the earth is rE; ϕ denotes the angle between the ray path (dashed line) and the radial direction r; and the bending angle, impact parameter, and tangent radius of a given ray path are represented by α, a, and rt, respectively.

Fig. 2.

Illustration of the GPS RO geometry. The radius of the earth is rE; ϕ denotes the angle between the ray path (dashed line) and the radial direction r; and the bending angle, impact parameter, and tangent radius of a given ray path are represented by α, a, and rt, respectively.

Normally, by using the substitution x = n(r)r and a = n(rt)rt, (3) can be transformed to (4) when x(r) is a monotonic function:

 
formula

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

 
formula

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 h2 and h3, or r2 and r3 (tangent height: h = rrE). We separate the portion of the true refractive index profile below h3 into three intervals:

  • A, the mixing layer: surface to h1 in which h and n are referred to as hA and nA,

  • B, the shadow layer: h1 to h2 in which h and n are referred to as hB and nB,

  • C, the SR layer: h2 to h3 in which h and n are referred to as hC and nC.

Fig. 3.

Illustrations of x(h), n(r), α(x), and n(x) for a profile with an SR layer at the top of the boundary layer. Solid curves indicate the true profile and dashed curves indicate the Abel-retrieved profile.

Fig. 3.

Illustrations of x(h), n(r), α(x), and n(x) for a profile with an SR layer at the top of the boundary layer. Solid curves indicate the true profile and dashed curves indicate the Abel-retrieved profile.

As seen in Fig. 3a, within the shadow layer and the SR layer (h1h3), h(x) is a multivalued function. As the tangent points of rays descends toward the height of critical refraction (i.e., h3), 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 h1 and h3 are both equal to x1, and the SR and shadow layer collapse to a singularity (i.e., an infinite large bending angle) at x1 in the bending angle profile (Kursinski et al. 2000). The Abel-retrieved refractive index profile is indicated by dashed lines, that is, (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 < x1, we have

 
formula

where a = ñ(t)t = nA(rt)rt. 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 of the Abel retrieval refers to a height below r3, whereas rt of the true profile refers to a height below r1. The limits, t and rt, must correspond to the same impact parameter to produce the same α(a) profile. Since the true and the Abel-retrieved refractive index profiles in Fig. 3b are identical above r3 (or equivalently h3), the corresponding parts of the integrals in (6) above r3 are identical. Thus, we need only consider the parts of the integrals below r3, and instead of (6), we can write

 
formula

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

 
formula

Next, by taking advantage of the fact that x and h uniquely define n [i.e., x = nr = n(h + rE)], 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., hC(x) and hB(x)] can be approximated by two straight-line segments:

 
formula
 
formula

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)

 
formula

where z = x1x/x2x1.

Given x1, x2, h1, and h2, (10) relates the Abel-retrieved refractivity profile (x) to a profile below h1, that is, hA(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., hB(x) and hC(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 h1 travel through these layers, and the bending angle below x1 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 (x1) and the upper limit of the SR layer (h3) 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 h3 and x1 could become a challenging problem and will be discussed in section 5. For simplicity, in the following we assume h3 and x1 are known parameters obtained from the bending angle and refractivity retrievals.

Knowing h3 and x1, (10) has two degrees of freedom, h1 and x2. Based on the simple linear segment assumption for hB(x) and hC(x), there is one more degree of freedom left for solving hB(x) and hC(x), namely, the lower limit of the SR layer, h2. To emphasize the degrees of freedom we write (10), (9a), and (9b) as follows:

 
formula
 
formula
 
formula

Besides the independent variable, x, hB(x), and hC(x) depend on x2, h1, and also h2, whereas hA(x) does not depend on h2 but only depends on x2 and h1. 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 x2, h1, and h2. However, not all combinations of x2, h1, and h2 make sense physically. For example, it is obvious from the formulation of the problem that x2x1, and that h1h2h3. It is interesting to point out that the second term on the rhs of (10), or equivalently Gh1, x2(x) in (11a), actually is a measure of the error in the Abel-retrieved refractivity profile, when hA(x) represents the “truth.” It is easily verified that Gh1, x2(x) (x0xx1) is always nonpositive, and so (x) ≥ hA(x) and ñ(x) ≤ nA(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 (h3h1), (x2x1), and (x1x), 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 (h3h1). Also note that the function −Gh1, x2(x) increases monotonically as z decreases. Accordingly, the larger the refractivity gradient inside the SR (i.e., x2x1), 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 hA(x) is finite and continuous at the bottom of the shadow layer, that is,

 
formula

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 x0 from the retrieved bending angle profile can be used as a surface constraint. Thus, as a second constraint we have

 
formula

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 (x0). Possible alternatives to the constraint in (13) will be discussed in section 5.

Equation (12) indicates that hB(x) ≈ hA(x) = h1 − (h2h1)z2 for small z (x close to x1). Therefore, Taylor expansion of (10) yields the following equation valid for small z:

 
formula

Invoking the definition of z = x1x/x2x1, we see that to the first order, (x) − h3 behaves like a square root function in x near x1. It follows that

 
formula

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

Each set of h1 and x2 thus solved for uniquely defines the structure of hA(x) according to (10). Because of the straight-line parameterization of hB(x) and hC(x), and the continuous slope constraint (11), we can determine h2 simply by linear extrapolation of hA(x) until x = x2. Thus hB(x) and hC(x) are obtained via (9a) and (9b). For different values of hA(x0), 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 x0 would be the lowest ray path with the tangent point touching the earth’s surface. Thus, the solution satisfying hA(x0) = 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 hA(x0) = 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.

Fig. 4.

End-to-end simulation system in the presence of SR (geometrical optics).

Fig. 4.

End-to-end simulation system in the presence of SR (geometrical optics).

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 h3, and its corresponding impact parameter x1, 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., hA(x0) = 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.

Fig. 5.

The Abel retrieval (dashed), a subset of a continuum (family) of refractivity solutions (dotted), and the reconstructed solution (dash–dot) in terms of (a) N(h) and (c) h(x) are consistent with (d) one bending angle profile. (b) The percentage refractivity errors for the Abel and reconstructed profiles compared to the true profile (solid).

Fig. 5.

The Abel retrieval (dashed), a subset of a continuum (family) of refractivity solutions (dotted), and the reconstructed solution (dash–dot) in terms of (a) N(h) and (c) h(x) are consistent with (d) one bending angle profile. (b) The percentage refractivity errors for the Abel and reconstructed profiles compared to the true profile (solid).

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.

Fig. 6.

Reconstruction results for four different MBL structures over global oceans with each row demonstrating one simulation case. (left) The sounding profile (solid), the Abel retrieval (dashed), and the reconstructed solution (dotted) in terms of x [=n(r)r] vs height. (middle) Similar to left panel, but in terms of refractivity vs height. (right) The percentage errors of the Abel retrieval (dashed) and the reconstructed solution (dotted) compared to the sounding profile. The three dashed horizontal lines in each panel represent h1, h2, and h3 upward.

Fig. 6.

Reconstruction results for four different MBL structures over global oceans with each row demonstrating one simulation case. (left) The sounding profile (solid), the Abel retrieval (dashed), and the reconstructed solution (dotted) in terms of x [=n(r)r] vs height. (middle) Similar to left panel, but in terms of refractivity vs height. (right) The percentage errors of the Abel retrieval (dashed) and the reconstructed solution (dotted) compared to the sounding profile. The three dashed horizontal lines in each panel represent h1, h2, and h3 upward.

Fig. 7.

Similar to Fig. 6, but for two cases not following very well the parameterization used in the reconstruction method. The three dashed horizontal lines in each panel represent h1, h2, and h3 upward.

Fig. 7.

Similar to Fig. 6, but for two cases not following very well the parameterization used in the reconstruction method. The three dashed horizontal lines in each panel represent h1, h2, and h3 upward.

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 h2, 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., h1).

Fig. 8.

Plot of refractivity errors in 55 cases for both Abel and the reconstructed profiles, offset by the bottom height of the shadow layer, h1.

Fig. 8.

Plot of refractivity errors in 55 cases for both Abel and the reconstructed profiles, offset by the bottom height of the shadow layer, h1.

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 h2. However, the dropsonde and radiosonde profiles normally bear a rounding edge at h2. Thus, the reconstruction method tends to overestimate the refractivity below h2 and underestimate it inside the SR layer. Also note that, due to the large refractivity gradient inside the SR layer, a small error in h2 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 (x1) and the upper limit of the SR layer (h3) 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 (x1 or h3) 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

REFERENCES
Albrecht
,
B. A.
,
V.
Ramanathan
, and
B. A.
Boville
,
1986
:
The effects of cumulus moisture transports on the simulation of climate with a general circulation model.
J. Atmos. Sci.
,
43
,
2443
2462
.
Ao
,
C. O.
,
T. K.
Meehan
,
G. A.
Hajj
,
A. J.
Mannucci
, and
G.
Beyerle
,
2003
:
Lower troposphere refractivity bias in GPS occultation retrievals.
J. Geophys. Res.
,
108
.
4577, doi:10.1029/2002JD003216
.
Beljaars
,
A. C. M.
, and
P.
Viterbo
,
1998
:
Role of the boundary layer in a numerical weather prediction model.
Clear and Cloudy Boundary Layers, A. A. M. Holtslag and P. G. Duynkerke, Eds., Royal Netherlands Academy of Arts and Sciences, 287–304
.
Betts
,
A. K.
, and
W.
Ridgeway
,
1989
:
Climate equilibrium of the atmospheric convective boundary layer over a tropical ocean.
J. Atmos. Sci.
,
46
,
2621
2641
.
Beyerle
,
G.
,
M. E.
Gorbunov
, and
C. O.
Ao
,
2003
:
Simulation studies of GPS radio occultation measurements.
Radio Sci.
,
38
.
1084, doi:10.1029/2002RS002800
.
Born
,
M.
, and
E.
Wolf
,
1964
:
Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light.
Pergamon Press, 856 pp
.
Bretherton
,
C. S.
, and
Coauthors
,
2004
:
The EPIC 2001 stratocumulus study.
Bull. Amer. Meteor. Soc.
,
85
,
967
977
.
Clarke
,
R. H.
, and
R. R.
Brook
,
1979
:
The Koorin Expedition—Atmospheric Boundary Layer Data over Tropical Savannah Land.
Australian Government Publishing Service, 359 pp
.
Clarke
,
R. H.
,
A. J.
Dyer
,
R. R.
Brook
,
D. G.
Reid
, and
A. J.
Toup
,
1971
:
The Wangara Experiment: Boundary-layer data. Tech. Paper 19, Division of Meteorological Physics, CSIRO, 21 pp
.
Cosma-Averseng
,
S.
,
C.
Flamant
,
J.
Pelon
,
S. P.
Palm
, and
G. K.
Schwemmer
,
2003
:
The cloudy atmospheric boundary layer over the subtropical South Atlantic Ocean: Airborne-spaceborne lidar observations and numerical simulations.
J. Geophys. Res.
,
108
.
4220, doi:10.1029/2002JD002368
.
Curran
,
R. J.
,
1989
:
Satellite-borne lidar observations of the Earth: Requirements and anticipated capabilities.
Proc. IEEE
,
77
,
478
490
.
Eshleman
,
V. R.
,
1973
:
The radio occultation method for the study of planetary atmospheres.
Planet. Space Sci.
,
21
,
1521
1531
.
Fjeldbo
,
G.
, and
V. R.
Eshleman
,
1968
:
The atmosphere of Mars analyzed by integral inversion of the Mariner IV occultation data.
Planet. Space Sci.
,
16
,
1035
1059
.
Fjeldbo
,
G.
,
A. J.
Kliore
, and
V. R.
Eshleman
,
1971
:
The neutral atmosphere of Venus as studied with the Mariner V radio occultation experiment.
Astron. J.
,
76
,
123
140
.
Garratt
,
J. R.
,
1992
:
The Atmospheric Boundary Layer.
Cambridge University Press, 334 pp
.
Gorbunov
,
M. E.
,
2002
:
Canonical transform method for processing radio occultation data in the lower troposphere.
Radio Sci.
,
37
.
1076, doi:10.1029/2000RS002592
.
Hajj
,
G. A.
, and
Coauthors
,
2004
:
CHAMP and SAC-C atmospheric occultation results and intercomparisons.
J. Geophys. Res.
,
109
.
D06109, doi:10.1029/2003JD003909
.
Heckley
,
W. A.
,
1985
:
Systematic errors in the ECMWF operational forecasting model in tropical regions.
Quart. J. Roy. Meteor. Soc.
,
111
,
709
738
.
Hong
,
S. Y.
, and
H. L.
Pan
,
1996
:
Nonlocal boundary layer vertical diffusion in a medium-range forecast model.
Mon. Wea. Rev.
,
124
,
2322
2339
.
Izumi
,
Y.
,
1971
:
Kansas 1968 field program data report. Air Force Cambridge Research Laboratory Paper 379, AFCRL-72-0041, Bedford, MA, 79 pp
.
Izumi
,
Y.
, and
S. J.
Caughey
,
1976
:
Minnesota 1973 atmospheric boundary layer experimental data report. Air Force Cambridge Research Laboratory Rep. AFCRL-TR-76-0038, Bedford, MA, 28 pp
.
Jensen
,
A. S.
,
M. S.
Lohmann
,
H-H.
Benzon
, and
A. S.
Nielsen
,
2003
:
Full spectrum inversion of radio occultation signals.
Radio Sci.
,
38
.
1040, doi:10.1029/2002RS002763
.
Klein
,
S. A.
, and
D. L.
Hartmann
,
1993
:
The seasonal cycle of low stratiform clouds.
J. Climate
,
6
,
1587
1606
.
Kursinski
,
E. R.
,
G. A.
Hajj
,
J. T.
Schofield
,
R. P.
Linfield
, and
K. R.
Hardy
,
1997
:
Observing Earth’s atmosphere with radio occultation measurements using the Global Positioning System.
J. Geophys. Res.
,
102
,
23429
23466
.
Kursinski
,
E. R.
,
G. A.
Hajj
,
S. S.
Leroy
, and
B.
Herman
,
2000
:
The GPS radio occultation technique.
Terr. Atmos. Oceanic Sci.
,
11
,
53
114
.
Lettau
,
H. H.
, and
B.
Davidson
,
1957
:
Exploring the Atmosphere’s First Mile.
Pergamon Press, 578 pp
.
McCormick
,
M. P.
, and
Coauthors
,
1993
:
Scientific investigations planned for the Lidar In-space Technology Experiment (LITE).
Bull. Amer. Meteor. Soc.
,
74
,
205
214
.
Mechoso
,
C. R.
, and
Coauthors
,
1995
:
The seasonal cycle over the tropical Pacific in coupled ocean–atmosphere general circulation models.
Mon. Wea. Rev.
,
123
,
2825
2838
.
Rocken
,
C.
, and
Coauthors
,
1997
:
Analysis and validation of GPS/MET data in the neutral atmosphere.
J. Geophys. Res.
,
102
,
29849
29866
.
Rocken
,
C.
,
Y. H.
Kuo
,
W.
Schreiner
,
D.
Hunt
,
S.
Sokolovskiy
, and
C.
McCormick
,
2000
:
COSMIC system description.
Terr. Atmos. Oceanic Sci.
,
11
,
21
52
.
Smith
,
E. K.
, and
S.
Weintraub
,
1953
:
The constants in the equation for atmospheric refractive index at radio frequencies.
Proc. IRE
,
41
,
1035
1037
.
Sokolovskiy
,
S. V.
,
2001a
:
Modeling and inverting radio occultation signals in the moist troposphere.
Radio Sci.
,
36
,
441
458
.
Sokolovskiy
,
S. V.
,
2001b
:
Tracking tropospheric radio occultation signals from low Earth orbit.
Radio Sci.
,
36
,
483
498
.
Sokolovskiy
,
S. V.
,
2003
:
Effect of superrefraction on inversions of radio occultation signals in the lower troposphere.
Radio Sci.
,
38
.
1058, doi:10.1029/2002RS002728
.
Swinbank
,
W. C.
,
1968
:
A comparison between predictions of dimensional analysis for the constant flux layer and observations in unstable conditions.
Quart. J. Roy. Meteor. Soc.
,
94
,
460
467
.
von Engeln
,
A.
, and
J.
Teixeira
,
2004
:
A ducting climatology derived from the European Centre for Medium-Range Weather Forecasts global analysis fields.
J. Geophys. Res.
,
109
.
D18104, doi:10.1029/2003JD004380
.
Winker
,
D.
,
R.
Couch
, and
M.
McCormick
,
1996
:
An overview of LITE: NASA’s Lidar In-space Technology Experiment.
Proc. IEEE
,
84
,
164
180
.
Zeng
,
X.
,
M. A.
Brunke
,
M.
Zhou
,
C.
Fairall
,
N. A.
Bond
, and
D. H.
Lenschow
,
2004
:
Marine atmospheric boundary layer height over the eastern Pacific: Data analysis and model evaluation.
J. Climate
,
17
,
4159
4170
.

APPENDIX

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

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

 
formula

Since n(x) = x/[rE + h(x)] and ñ(x) = x/[rE + (x)], and because hrE, (A1) can be approximated as

 
formula

Using y instead of a, and writing Δh(x) = hC(x) − hB(x), we define the function f (y) as

 
formula

Thus, from (A2) and (A3) we have

 
formula

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

 
formula

Having hA(x1) = h1 and (x1) = h3 (cf. Fig. 3a), we arrive at

 
formula

Equation (A4) is a general expression relating hA(x′) to (x′), given h1, h3, 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

 
formula

Integrating the last expression by parts yields

 
formula

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 = hChB = A(x2x), where A = (h3h1)(x2x1)−1. We may now write (A5) as

 
formula

With t = x2x21, the second term on the rhs of (A6) can be written as

 
formula

Using the substitution t = tan−1(x2x21/x21x2), the third term on the rhs of (A6) can be written as

 
formula

Since x21x2x21, 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

 
formula

Substituting (A7) and (A9) into (A6), and then substituting (A6) into (A4), we arrive at

 
formula

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

 
formula

where z = x1x/x2x1.

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