Interactive comment on “ Ionospheric correction of GPS radio occultation data in the troposphere ” by Z

This paper considers the ionospheric correction of GPS radio occultation in the troposphere, and how to extrapolate the (alpha_l1-alpha_l2) differences. The information in the paper will be of interest to other GPS-RO data providers, and NWP users who may wish assimilate the data. However, I recommend publish after major revision, because often the results are explained on the basis of simulations that are not presented in any detail. For example, the "ionosphere induced vertical shifts" noted in the abstract, are similar to shifts seen in wave optics modeling. But this modeling is only mentioned, and not presented in any meaningful detail. Hence, it is difficult for the reader to understand how the electron density at the LEO can produce this result (particularly when the ray tracing results differ).


Introduction
When GPS radio occultation (RO) signals are used for monitoring the neutral atmosphere, the ionospheric effect has to be removed.While the neutral atmospheric effect on the RO signals exponentially decreases with the height, the ionospheric effect on the average slowly increases with the height.Thus, removal of the ionospheric effect (i.e., ionospheric correction) is most important in the stratosphere and above.Nevertheless, neglecting the ionospheric correction in the troposphere results in a small but statistically significant inversion bias.One of the common methods of dual-frequency, modelindependent ionospheric correction is a linear combination of the L1 and L2 GPS RO observables and, in particular, the linear combination of the L1 and L2 bending angles taken at the same impact parameter (Vorob'ev and Krasil'nikova, 1994) (hereafter the "standard ionospheric correction").However, a disadvantage of the standard ionospheric correction is the amplification of uncorrelated noises (errors) on L1 and L2.
When the first GPS RO data became available in 1995 (Ware et al., 1996), it became clear that the standard linear combination of the L1 and L2 bending angles is useless in the troposphere, mainly due to noise and tracking errors on the encrypted L2P signal.In order to approximately remove the mean ionospheric effect, it was proposed that the L1 bending angle be corrected by using the L1-L2 bending angle extrapolated from above the troposphere (Kursinski et al., 1997(Kursinski et al., , 2000;;Rocken et al., 1997;Steiner et al., 1999;Hajj et al., 2002).As a result of the open-loop (OL) tracking of the L1 in the troposphere (Sokolovskiy et al., 2009b;Ao et al., 2009), the L2P is unavailable and only the ionospheric correction by extrapolation of the L1-L2 can be applied.With the OL tracking of the un-encrypted L2C (Sokolovskiy et al., 2014), both the L1CA and L2C signals, free of the tracking errors, are available in the troposphere.Theoretically, this may allow an extension of the standard ionospheric correction into the troposphere.However, in practice, as noted by Steiner et al. (1999), not only the noise but also the small-scale refractivity irregularities in the troposphere may introduce additional errors into the standard ionospheric correction.As follows from the analysis of the occultations with the L1CA and L2C available in the OL mode (examples shown below in Sect.1), the standard ionospheric correction down to the surface is feasible for only smooth RO signals under lowmoisture conditions (e.g., high latitudes, local winter).In the presence of moisture, especially in the tropics, fluctuations Published by Copernicus Publications on behalf of the European Geosciences Union.
of the RO signals result in an increase of noise (errors) after the standard ionospheric correction.Also, as shown in this study, the horizontal inhomogeneity of the ionosphere may result in additional errors of the standard ionospheric correction around strong inversion layers.The ionospheric correction by extrapolation of the L1-L2 reduces the errors of these types.
For an optimal application of the ionospheric correction in the troposphere, it is necessary to define the transition height for replacement of the standard linear combination by extrapolation of the L1-L2.It is also necessary to define the extrapolation function and the height interval in which this function is fitted to the L1-L2 and then used to model the L1-L2 below the transition height.In previous studies (Rocken et al., 1997;Steiner et al., 1999;Hajj et al., 2002), the transition height and the fitting interval were not always defined and simple extrapolation functions (constant and linear) were used primarily.In this study, the extrapolation function, which models the effects of the ionospheric F and E layers, is introduced, and a special case of this function is used for processing and statistical analysis of the GPS RO data.The selection of the transition height is based on the noisiness of the ionosphere-free bending angle.A transition too high results in an increase of noise due to the uncorrected, small-scale ionospheric effects.A transition too low results in an increase of noise due to the small-scale tropospheric irregularities and L2P tracking errors.An optimal transition height minimizes the noise of the ionosphere-free bending angle at all heights.As follows from the results of this study, the optimal transition height depends on the latitude and is different for different GPS RO receivers and firmware, for L2P and L2C signals, and for setting and rising occultations.
Section 2 presents examples of the occultations with L1CA and L2C available in the OL mode and with the standard ionospheric correction applied down to the surface.Section 3 introduces a model of the L1-L2 bending angle for extrapolation of the ionospheric correction in the troposphere.Section 4 presents statistical distributions of the dynamic transition heights (as determined individually for each occultation).Section 5 presents the statistical comparison of the GPS RO bending angles to those from the European Centre for Medium-range Weather Forecasting (ECMWF) forecast for different static transition heights.This allows us to determine optimal static transition heights.Section 6 concludes the study.

Examples of the standard ionospheric correction in the troposphere
The effect of uncorrelated L1 and L2 random errors on the ionospheric correction is known.Let α 1 and α 2 be the L1 and L2 bending angles and β 1 and β 2 their observational noises (errors), whereas f 1 = 1.57542GHz and f 2 = 1.2276GHz are the GPS frequencies.The 1st-order ionosphere-free bend-ing angle is equal to where h = a − r c − h g is the impact height (a is the impact parameter, r c is the local curvature radius of reference ellipsoid, and h g is the geoid undulation), For the uncorrelated, random errors β 1 and β 2 , the root mean square (rms) error of the ionosphere-free bending angle equals If the uncorrelated L1 and L2 random errors have the same rms magnitudes (< β 2 1 >= < β 2 2 >), then, after the ionospheric correction, the rms error is amplified by a factor of (c 2 1 + c 2 2 ) 1/2 ∼ = 3.If the L2 random error is much larger than the L1 error, then it is feasible to correct α 1 by an additionally smoothed α 1 − α 2 : (3) This reduces the effect of L2 errors but increases the errors due to un-corrected small-scale ionospheric effects on L1.The optimal smoothing interval for α 1 − α 2 is determined by minimizing the combined effect of both errors individually for each occultation (Sokolovskiy et al., 2009a).This approach is routinely used in the COSMIC Data Analysis and Archive Center (CDAAC) processing (Schreiner et al., 2011).
Figure 1 shows α 1 (h), α 2 (h), and α(h) for four Constellation Observation System for Meteorology, Ionosphere, and Climate (COSMIC) occultations with the L1CA and L2C acquired in the OL mode in the troposphere.Bending angles were calculated from the complex RO signals with the use of a wave optics (WO) transform (phase matching) (Jensen et al., 2004) and then smoothed with a window of 0.1 km. Figure 1a shows a high-latitude wintertime occultation.Both α 1 (h) and α 2 (h) are smooth functions due to low humidity resulting in small fluctuations of refractivity.The fluctuation of α(h) is not substantially different from those of α 1 (h) and α 2 (h).The standard ionospheric correction performs satisfactorily in the troposphere.
Figure 1b shows a tropical occultation (from the center of the Pacific Ocean).Commonly, most of the tropical occultations are affected by strong random refractivity irregularities caused by moist convection in the troposphere, which broaden the spectra of RO signals both in time and impact height representations (Gorbunov et al., 2006;Sokolovskiy et al., 2010).Correspondingly, α 1 (h) and α 2 (h) have strong fluctuations that are partially uncorrelated (Sokolovskiy et al., 2014).After applying the standard ionospheric correction, these fluctuations are amplified as seen in Fig. 1b.Thus, in the moist convective troposphere, the standard ionospheric correction, generally, results in substantial increase of noise.Such increased errors can be often (but not always) observed around the inversion layers in L2C occultations.They can be viewed as resulting from different shifts of the impact heights at L1 and L2 frequencies.An explanation and modeling of these errors is discussed in Appendix A. Extrapolation of the ionospheric correction (discussed in the next section) reduces the errors such as those below ∼ 7 km in Fig. 1b and below ∼ 3 and ∼ 3.5 km in Fig. 1c and d.

Ionospheric correction by extrapolation of L1-L2 bending angle
At heights where errors of α 2 (h) become too large (such as in the troposphere), it makes sense to approximate α 1 (h) − α 2 (h) by a smooth function α ext (h), which is fitted to α 1 (h)− α 2 (h) at some interval (h ext , h top ) where the errors of α 2 are smaller: and, at h < h ext , replace Eq. (3) with The ionospheric correction by extrapolation eliminates the errors of α 2 by instead introducing the errors due to uncorrected small-scale ionospheric effects on α 1 .If the α 2 errors are larger than the small-scale ionospheric effects, the ionospheric correction by extrapolation results in the overall reduction of the magnitude of bending angle error.However, such replacement of the errors may also change the vertical error correlation, as noted by Syndergaard et al. (2013).Though detailed investigation of this effect is outside the scope of this study, it is important to note that increase of the error correlation radius increases the error of retrieved refractivity for a given rms bending angle error (see, for example, frequency response of the Abel inversion in Lohmann (2005), Fig. 1).Thus reduction of the rms magnitude of the bending angle error, by itself, does not warrant reduction of the refractivity error.Since, on average, α 1 − α 2 is a smooth function, in previous studies, α ext was either a constant or linear function (Kursinski et al., 1997(Kursinski et al., , 2000;;Rocken et al., 1997;Steiner et al., 1999;Schreiner et al., 2011) fitted to the observational α 1 − α 2 in some interval above the extrapolation height.However, α 1 − α 2 is also affected by the ionospheric irregularities and L2 tracking errors.To reduce these effects on α ext , it makes sense to increase the fitting interval and to apply a more complicated model of the ionospheric effect on bending angle.For this purpose, we use an approximate expression for the bending angle response δα from an infinitely thin refractivity layer at a height z 0 : where δ(z) is a delta function.The ionospheric layer, such as F2, is not thin.Figure 2 shows electron density profiles modeled by Chapman layer (left panel) and α 1 (h) − α 2 (h) calculated by ray tracing and approximated by the Eq. ( 6) fitted in the interval 20 km < h < 80 km (right panel).Details are provided in the figure caption.It is seen that the model (Eq.6) satisfactorily represents α 1 (h) − α 2 (h) for modeling both F2 layer at height 300 km and E layer at height 100 km.
The accuracy of the approximation reduces with the increase of the thickness of the layer, which is to be expected.We model the ionospheric bending angle by the function where z E = 100 km and z F2 = 300 km are approximate heights of the E and F2 layers and search for the coefficients A, B, C, and D by least-squares fitting α ext (h) to the obser- and α ext (h) for six COSMIC occultations for h ext = 20 km. Figure 3b shows the same but with an excluded fourth term in Eq. ( 7).It is seen from comparison of Fig. 3a and b that keeping the fourth term in α ext (h) results in an overfitting.Our statistical analysis (details are omitted) also shows strong cross-correlation between the coefficients B and D. Thus it is sufficient to keep the third term, which models the response from E layer, while the response from F2 layer, as well as the effects of horizontal inhomogeneity of the ionosphere (see Appendix A), are modeled by the first two terms: We found that the use of the fitting function (Eq.8), compared to the fitting functions with only the first term (constant) or first two terms (linear fit) (fitted in a shorter interval, 10 km above the transition height), results in slightly smaller standard deviation of RO bending angles from those for the collocated ECMWF analyses.In the following, we use the function (Eq.8) for extrapolation of the ionospheric correction.

Dynamic transition height for extrapolation of the ionospheric correction
In one of the options used in the CDAAC inversion algorithm, the L2 quality is checked from top to bottom based on a simple criterion | s 1 − s 2 | < 6 cm, where s 1 and s 2 are raw (unsmoothed) excess phase changes between 50 Hz samples (Kuo et al., 2004).The standard ionospheric correction is replaced by extrapolation below the height (called the L2 drop height, same as the transition height), defined as the Figure 4 shows the latitudinal distributions of the L2 drop heights for April 2012.Figure 4a and b show the L2P drop heights for the COSMIC setting and rising occultations.The upper part of the "donut-shaped" structure in Fig. 4b is correlated with the height of the tropopause.In many occultations, the sharp structure of the tropopause is sufficient to cause a strong enough fluctuation of the RO signal resulting in a tracking instability of the L2P.For those setting occultations where the L2P tracking remains sufficiently stable at the tropopause, it becomes unstable in the presence of fluctuations caused by the tropospheric moisture, or the L2 becomes unavailable below the transition height from the phase-locked loop (PLL) to the OL mode; this explains the lower part of the donut-shaped structure.For rising COS-MIC occultations, the structure of the distribution of the L2P drop heights is, in general, similar to that of the setting occultations, except that the lower part is higher because the OL-PLL transition generally needs extra time for locking on the L2P and this, on average, happens higher than the PLL-OL transition for setting occultations.The L2C drop heights for COSMIC setting occultations (Fig. 4c) are primarily related to the PLL-OL transition; this confirms that the L2C PLL tracking is stable and not susceptible to signal fluctuations (Sokolovskiy et al., 2014).Although there are some differences, the structures of the distributions of the L2P drop heights for the Meteorological Operational satellite program (Metop) occultations prior to the 2013 firmware update (Figs.4d and e) are, in general, similar to those for COSMIC.For setting occultations (Fig. 4e), fewer occultations are affected by the L2P tracking instability induced by the tropopause, and, due to the dynamic PLL-OL transition, the L2P can, on average, be more stably tracked down to the lower height in the troposphere than for COSMIC (Fig. 4b).For rising occultations (Fig. 4d), the lock on the L2P happens, on average, at lower heights than for COS-MIC (Fig. 3a).Overall, the comparison of Fig. 4a and b to d and e suggests less noisy (more stable) L2P tracking for Metop than for COSMIC.
As follows from Fig. 4, the dynamic approach for extrapolation of the ionospheric correction results in a rather large spread of the L2P drop heights in low latitudes as well as systematic differences of the L2 drop heights for different missions, rising and setting occultations, and L2P and L2C signals.Generally, this approach may be optimal for the processing of RO data for general purposes and weather applications.An alternative (static) approach applies extrapolation of the ionospheric correction below a fixed transition height that is predetermined based on statistical evaluation and minimization of the inversion errors.We believe that the static approach, when all occultations are processed in the same way, may be superior for climate applications.

Estimation of the optimal fixed transition height for extrapolation of the ionospheric correction
In this section, we investigate the effect of different fixed transition heights for extrapolation of the ionospheric correction on the retrieved ionosphere-free bending angle profiles.
For retrieval of the L1 bending angle, the phase matching method (Jensen et al., 2004) is applied up to 20 km or the transition height h ext (whichever is higher), and geometric optics is applied above.For retrieval of the L2 bending angle, the geometric optics is applied.To evaluate the retrieved ionosphere-free bending angle profiles, they are statistically compared to the bending angles obtained from the ECMWF global forecasts, by calculating the mean (bias) and standard deviation.A transition too high results in the increase of random errors due to uncorrected, small-scale ionospheric effects, while a transition too low results in the increase of random errors due to noise and tracking errors on the L2.Although the ECMWF model, which is used as the reference, also has errors, it is important that both the errors related to uncorrected ionospheric effects and L2 noise are uncorrelated with the model forecast errors.Consequently, the minimal standard deviation (where all independent errors are summed with squares) can be used as the criterion for finding an optimal transition height.The mean deviation (where all biases are summed with their signs), generally, cannot be used as such a criterion unless the bias becomes too large and can be clearly attributed to RO rather than to the model.transition heights applied for extrapolation of the ionospheric correction.In each case, an optimal transition height is located between those fixed heights which result in smaller standard deviations at all heights.More accurately, the optimal transition height can be estimated as the height where the difference between these standard deviations for the fixed heights changes the sign.
Figure 5 shows the statistics for the COSMIC L2P occultations with transition heights at 25, 20, 15, and 10 km over three latitude bands.It is apparent that if the transition is too high, such as at 25 km, the standard deviation in some interval (∼ 5 km) below the transition height increases, which indicates the noise from uncorrected small-scale ionospheric structures.If the transition is too low, such as at 15 or 10 km, the standard deviation in some interval above the transition height (∼ 5-10 km) substantially increases.Clearly, involving the noisy L2P in the ionospheric correction significantly deteriorates the statistical results.The optimal transition heights for different latitudes are not very different, ranging from ∼ 17 km at high latitudes to ∼ 19 km at low latitudes.Though the vertical error correlation might change due to the ionospheric correction by extrapolation, the statistical results for refractivity (right panels) affected by error propagation through the Abel transform, are overall consistent with those for bending angle (left panels).Figure 6 shows the statistics for the COSMIC L2P rising and setting occultations.Again, the optimal transition heights are not very different: slightly lower for setting (∼ 17 km) than for rising (∼ 19 km) occultations.This difference may suggest, on average, lower quality of the L2P signal right after lock for rising occultations than after tracking for an extended time for setting occultations.We also examine the statistics for the Metop L2P rising and setting occultations before the 2013 firmware update in Fig. 7.The optimal transition heights are ∼ 15 and ∼ 20 km for setting and rising occultations, respectively.A comparison of Figs. 6 and 7 shows that Metop has a smaller standard deviation (i.e., better L2P quality) than COSMIC, especially for setting occultations.Figure 8a and b show statistics for the COSMIC L2C setting occultations (rising L2C occultations currently are not available).The optimal transition height is about 10 km.Because the L2C is tracked in the OL mode down to the bottom of the occultations, Fig. 8c and d also show the standard deviation in the lower troposphere for the transition height at 10 km and for the standard ionospheric correction without extrapolation.It is apparent that application of the standard ionospheric correction down to the bottom of the occultations makes the ionosphere-free bending angles in the lower troposphere, on average, noisier than the correction by extrapolation below 10 km.

Conclusions
Application of the standard ionospheric correction for the L1 and L2 GPS RO bending angles results in an increase of random errors in the troposphere due to noise and tracking errors (mainly on the L2P signal) and tropospheric irregularities, and is also limited by availability of the L2 signal.Correction of the L1 bending angle by the L1-L2 bending angle extrapolated from above reduces these errors by introducing other errors due to uncorrected small-scale ionospheric effects.The extrapolation height, which minimizes the combination of both errors, is considered optimal and depends mainly on the of the L2 signal, which, in turn, depends on the receiver, tracking firmware, signal strength, and structure, and may depend on the tropospheric and ionospheric irregularities.
Results of this study show that the optimal transition height for the COSMIC L2P occultations is about 20 km (the height differs slightly for different latitudes and for rising and setting occultations).For Metop occultations before the 2013 firmware update, the optimal transition height is about 20 km for rising and ∼ 15 km for setting occultations, indicating a less noisy L2P than for COSMIC.
For the COSMIC L2C setting occultations (tracked in the OL mode in the lower troposphere), the optimal transition height is about 10 km.Extension of the standard ionospheric correction without extrapolation down to the bottom of the occultations, generally, increases the noisiness of the ionosphere-free bending angles in the troposphere.In particular, most significant increase is observed (i) in the presence of non-spherically symmetric refractivity irregularities (such as the moist convection) and (ii), sometimes, around strong inversion layers (such as the top of boundary layer).The latter effect is modeled and explained by horizontal in- homogeneity of the ionosphere.The results of the modeling suggest the ionospheric correction of the boundary layer depth determined from RO bending angles (application of this correction is a subject of a separate study).To minimize the noise of the ionosphere-free bending angles, occultations from different missions, rising and setting occultations, and L2P and L2C occultations can all be processed with different extrapolation heights.Alternatively, an optimal extrapolation height can be determined for each occultation individually.This may be better for weather applications.However, when GPS RO data are used for climate applications, different or dynamically determined extrapolation heights may result in different biases, which may propagate into climate signals.Thus, for climate applications, a constant extrapolation height may be superior.Based on the results of this study, for currently available GPS RO data, a transition height of 20 km (when the L2 is available down to that height) may be considered reasonable.

Appendix A: Modeling of the errors of the ionospheric correction induced by horizontal gradients in the ionosphere
The ionospheric correction Eq. ( 1) is based on the assumption that the observed α(h) can be decomposed as follows: where α atm (h) and α ion (h)/f 2 are to the bending angle from the neutral atmosphere and the ionosphere, respectively.This is true for a spherically symmetric refractivity and non-overlapping atmosphere and ionosphere because the bending angle can be represented in the form of an integral and split into two parts (Vorob'ev and Krasil'nikova, 1994;hereafter VK94).It is important that in this case α 1 (h) − α 2 (h) only depends on the ionospheric refractivity.
In case of non-spherically symmetric refractivity, the decomposition Eq. ( A1) may not always be applicable.Thus the standard ionospheric correction Eq. ( 1), generally, may result in errors.In this case, α and h derived from RO Doppler, generally, are different from the true bending angle and impact height (Healy, 2001).The α 1 (h) − α 2 (h), as well as the errors of the ionospheric correction, may depend on the atmospheric structure.In this section, by applying ray tracing and realistic modeling of α(h) derived from Doppler, we investigate the increase of the errors of ionospheric correction around strong inversion layers in the troposphere in the case of non-spherically symmetric ionosphere.We describe our modeling with the level of details sufficient for reproducing the results.Comparison with the results VK94 is discussed at the end of the section.A layout of ray tracing is shown in Fig. A1.A transmitter (GPS) and receiver (in low Earth orbit: LEO) are in co-planar circular orbits.This simplifies calculations and is sufficient to model and explain the errors of the ionospheric correction.The Earth radius is r e = 6370 km, and GPS and LEO orbit radii are r 1 = 26 600 km and r 2 = r e + 730 km = 7100 km, GPS and LEO velocities are v 1 = 4 km s −1 and v 2 = 8 km s −1 .For the ray tracing, we use the ray equation in the form (Kravtsov and Orlov, 1990) where r is the radius vector, dτ = dσ/n, dσ is the differential of length, and n is the refractive index.We integrate the Eq.(A2) in Cartesian coordinates (x, y) while specifying the refractivity in polar coordinates (r, θ ) by applying conversion at each integration step.For the integration, we apply the Runge-Kutta method of the 4th order with the integration step ranging from 1 to 10 km.Above the ionosphere (r > 7870 km in our modeling), the ray is replaced by a straight line.We start each ray from r 1 at a fixed central angle and at a different zenith angle φ 1 and integrate until r = r 2 after passing the minimum r min .This approximately reproduces a realistic GPS-LEO observation geometry (because GPS is moving slower and is located at a larger distance than LEO).The last integration step requires iterations because LEO is in the ionosphere.We note that the iterations are not required to integrate the rays from LEO to GPS (because GPS is in a vacuum), but then, to reproduce a realistic GPS-LEO observation geometry, each ray must be started from a different central angle.At the end of each ray, we calculate the zenith angle φ 2 from the ray direction at r 2 .We also store the central angle θ 12 .We note that φ 1 and φ 2 are true zenith angles of the ray which are not known from RO observation.A RO observable is the frequency f d of the received signal with the carrier frequency f , where where n 1 and n 2 are refractive indices at GPS and LEO and c is the light velocity in a vacuum (relativistic terms are neglected in Eq.A3).Next we solve Eq. (A3) for the ray zenith angles, by assuming n 1 = n 2 = 1 and spherically symmetric refractivity, thus adding Snell's equation r 1 sin φ 1 = r 2 sin φ 2 .This yields We note that φ 1 and φ 2 , generally, are different from the true φ 1 and φ 2 .The bending angle α and impact parameter a are Thus obtained α and a, generally, are not equal to the true values of α and a. Healy (2001) investigated errors of bending angle and impact parameter induced by horizontal gradients.However, Eqs. (A4) and (A5) can be viewed as the conversion of the observable f d (t) (where t is time), with the use of the orbit data, into a new observable α( a).Once this observable α( a) is adequately modeled (i.e., similarly to the processing of RO observational Doppler data, as discussed above), the fact that α = α and a = a does not matter and is not an error source.However, an adequate modeling is not always possible.In particular, in our study we investigate the effect of the ionosphere on the errors of the ionospheric correction.In practice, modeling the ionospheric state is more difficult than of the neutral atmospheric state and is not commonly applied.Thus in our modeling we attribute the ionospheric effect on α( a) to the error.In the description of our modeling below we follow calculations (Eqs.A4 and A5) and, at the end, omit ∼.
We calculate a set of rays with the increment φ 1 = 2×10 −7 rad (this provides a sufficiently dense grid of impact heights in the lower troposphere).We calculate observation functions α 1 (a) and α 2 (a) independently for f 1 and f 2 GPS frequencies.Since α 1 (a) and α 2 (a) are specified on different grids, we interpolate them (by cubic spline) on the same impact parameter grid for the purpose of the ionospheric correction.
In our modeling of refractivity we use the following function (continuous with derivatives): w(z) = w(z − z 0 , z), where The refractivity model is composed of the atmospheric and ionospheric terms: The atmospheric model is the exponential function with the inversion layer: where N 0 = 300, H = 7 km, z 0 = 1.5 km, z = 0.1 km, and c = 0.05.The N atm (z) is shown in Fig. A2a.The N atm (z) results in the gradient of about −100 km −1 at z = 1.5 km, which is below critical (−154 km −1 ) but still results in a rather sharp maximum bending angle of about 0.03 rad (similar to those observed in Figs.1c and d).
The ionospheric model is decomposed into radial and horizontal factors as follows: where N e = 10 6 cm −3 is the norming electron density.Function R(z) is where z 1 = 200 km, z 1 = 100 km, z 2 = 450 km, and z 2 = 150 km.This results in the height of the maximum First, we calculate α atm (h) for only the atmospheric (spherically symmetric) refractivity N atm (z) and use it as the reference.Next, we calculate α 1 (h), α 2 (h), and α(h) = c 1 α 1 (h) − c 2 α 2 (h) for the atmospheric and the ionospheric spherically symmetric refractivities by setting T (θ ) = 1.The results are shown in Fig. A3a.From the zoomed section of Fig. A3a, it is seen that the ionosphere results in the effects on α 1 (h) and α 2 (h) which are eliminated in α(h) (α(h) is indiscernible from α atm (h)), which is to be expected.
Next we model the effect of horizontal inhomogeneity of electron density in the ionosphere.Figure A4 (lower panel) shows a set of traced rays, in the coordinates (r −r e , θ ) spanning the troposphere (θ is counted from the direction to GPS).Upper panel shows functions T (θ ) resulting in change of maximum electron density from 5×10 5 to 1.5×10 6 cm −3 over a horizontal distance of about 1300 km.Case B corresponds to θ 0 = 1.05 rad (rays at the height H max on the side of GPS).The bending angles are shown in Fig. A3b.From the zoomed section of Fig. A3b, it is seen that the ionospheric effects on α 1 (h) and α 2 (h) are eliminated in α(h) above the inversion layer by leaving some residual error below.Cases C and D correspond to θ 0 = 1.65 rad (rays at the height H max on the side of LEO) but have different signs ±(θ − θ 0 ).In these cases, the errors of the ionospheric correction are very large at the inversion layer.Furthermore they have different structures related to different signs of the horizontal refractivity gradient in the ionosphere.One more case E corresponds to θ 0 = 1.35 rad (rays in the troposphere).In this case, though the ionosphere is horizontally inhomogeneous, the results of the ionospheric correction of L1 and L2 bending angles are indiscernible from the spherically symmetric case A and thus are not shown in Fig. A3.It is seen that the effect of the horizontal ionospheric inhomogeneity on the errors of the ionospheric correction is substantially larger when the inhomo- geneity is located on the side of LEO rather than on the side of GPS.This can be explained by the geometry (larger distance to GPS and its slower motion compared to LEO) and by differential bending of rays in the troposphere.As seen from lower panel of Fig. A4, the cross section of the set of rays spanning the troposphere is much larger on the side of LEO than on the side of GPS.This difference should be responsible for the difference in magnitudes of the errors of the ionospheric correction.For different observational geometry, the errors may be different.
To further justify the statements made at the beginning of this section, Fig. A5 shows α 1 (h) − α 2 (h) and the errors of the ionospheric correction α 1 (h)−α 2 (h).Panels a, c1, and d1 show α 1 (h)−α 2 (h) for the atmospheric refractivity with (red lines) and without (c = 0) (black lines) the inversion layer.For the spherically symmetric ionosphere (panel a) there is no difference in α 1 (h) − α 2 (h) with and without the inversion layer; there is no error of the ionospheric correction (also follows from Fig. A3a).The noise on α 1 (h) − α 2 (h) in the presence of the inversion layer is related to numerical noise of ray tracing, which increases in the regions of strong refractivity gradients.The magnitude and slope of α 1 (h) − α 2 (h) are consistent with those modeled in Fig. 2. Panels c1 and d1 show α 1 (h) − α 2 (h) for the models of the ionospheric inhomogeneity C and D (shown in Fig. A4).In this case, there are large differences in α 1 (h) − α 2 (h) with (red lines) and without (black lines) the inversion layer.Furthermore, the struc- tures of α 1 (h) − α 2 (h) without the inversion layer are quite different from those for a spherically symmetric ionosphere: the absolute magnitude decreases upward, and the sign can be both negative and positive.This provides an additional justification for keeping linear terms in the fitting function α ext (h) (see Sect. 3).Panels c2 and d2 show errors of the ionospheric correction α 1 (h) − α 2 (h) for the inversion layer and the models of the ionospheric inhomogeneity C and D. Red lines show the errors of the standard ionospheric correction.Green lines show the errors of the correction by extrapolation by using α 1 (h) − α 2 (h) without the inversion layer (black lines in panels c1 and d1) as the α ext (h).It is seen that extrapolation of the ionospheric correction results in some reduction of the errors.
Alternatively, the structure of the ionospheric effects on α 1 (h) and α 2 (h) around strong inversion layers induced by horizontally inhomogeneous ionosphere can be interpreted as impact height shifts h 1 and h 2 (as opposed to the additive effects which are eliminated by the standard ionospheric correction).These shifts have different signs related to different signs of the horizontal refractivity gradient in the ionosphere.Such an interpretation may have a useful practical implication.The height of a large lapse of the L1 RO bending angle has been used as a tool to monitor the depth of the atmospheric boundary layer and its variability (Sokolovskiy et al., 2011;Ho et al., 2015).Since the impact height shifts induced by horizontally inhomogeneous ionosphere may be quite large, of about 100 m and more, they will affect the estimated boundary layer depth.In particular, this effect may be most significant for evaluation of diurnal variability of the boundary layer depth because the ionosphere has a strong diurnal cycle which, potentially, may be aliased into the atmospheric diurnal cycle.While the shifts h 1 and h 2 cannot be measured directly, their difference h 1 − h 2 can be mea- sured (as can be seen from Figs. 1c and d).This can be done at a given α corresponding to maximum gradient |dα/dh|, or averaged over some region of large gradient.Our modeling (Fig. A3c and d) shows that with a sufficient accuracy h 1,2 ∼ 1/f 2 1,2 .This allows a simple correction for h 1,2 based on the measured difference: h 1,2 = c 2,1 ( h 1 − h 2 ).Naturally, such a correction is possible for only L2C occultations, and its practical implementation requires additional study.
The effect of the horizontal inhomogeneity of the ionosphere on the accuracy of the ionospheric correction was also investigated in VK94.By applying ray tracing, it was concluded that the error is negligible.The difference with our conclusions can be explained by two reasons.First, VK94 used the model of horizontal ionospheric inhomogeneity which is similar to our case E in Fig. A4 for which we also did not find errors.Second, and most important, is that VK94 applied the ionospheric correction for the true bending angle as a function of the height of the ray tangent point.According to VK94, this observable can be decomposed into the atmospheric and ionospheric terms, similarly to Eq. (A1), even when the ionospheric refractivity is non-spherically symmetric.We verified the results of VK94 with ray tracing and confirmed that for this observable the errors of the ionospheric correction are negligible not only in case E but also in cases C and D. But this observable is not realistic because it cannot be derived from RO Doppler in a general case of non-spherically symmetric refractivity.For the realistic RO observable given by Eqs.(A4, A5), the errors of the ionospheric correction can be significant, as demonstrated in this section.

Figure
Figure 1c and d show two occultations affected by strong refractivity gradients on top of the boundary layer, resulting in the large bending angle gradients and lapses around h ∼ 2.8 km and h ∼ 3.3 km.In these regions, the errors of the standard ionospheric correction are clearly increased.Such increased errors can be often (but not always) observed around the inversion layers in L2C occultations.They can be viewed as resulting from different shifts of the impact heights at L1 and L2 frequencies.An explanation and modeling of these errors is discussed in Appendix A. Extrapolation of the ionospheric correction (discussed in the next section) reduces the errors such as those below ∼ 7 km in Fig.1b and below

Figure 2 .
Figure 2. Left panel: examples of electron density profiles modeled by Chapman layer.Blue lines: F2 layer at 300 km with vertical scale 50 km (dashed line) and 75 km (solid line).Green lines: E layer at 100 km with vertical scale 5 km (dashed line) and 7.5 km (solid line).Right panel: L1-L2 bending angles obtained by ray tracing for the profiles in left panel (blue and green lines); approximations by bending angles calculated for refractivities modeled by delta functions at 300 and 100 km (red lines).Numbers (%) show fractional rms deviations (averaged below 20 km) between L1-L2 bending angles and their approximations.

Figure 3 .
Figure 3. Examples of COSMIC L1-L2 bending angle profiles (thin lines) and fitting functions (thick lines).(a, b) correspond to different models of the fitting function.

Figure 4 .
Figure 4. Latitudinal distributions of the dynamic L2 drop heights for COSMIC and Metop occultations.
Figures 5 through 8, in left panels, show statistical comparisons of the RO-retrieved ionosphere-free bending angles to the bending angles forward-modeled from the ECMWF forecasts for April 2012.Because extrapolation of the ionospheric correction may change the bending angle error correlation which affects the refractivity error (mentioned in Sect.3), we also perform statistical comparison of refractivities.These results are shown in right panels of Figs.5-8.Solid and dashed lines show the mean and standard deviations, respectively.Different colors correspond to different

Figure 5 .
Figure 5. Mean differences (solid) and standard deviations (dashed) of retrieved COSMIC L2P bending angles (left) and refractivities (right) for different extrapolation (transition) heights (colors indicated in panel f), relative to the colocated ECMWF forecast data over low (a, b), mid-(c, d), and high (e, f) latitudes.

Figure 6 .
Figure 6.Similar to Fig. 5, except for the COSMIC L2P rising (a, b) and setting (c, d) occultations, for all latitudes.

Figure 7 .
Figure 7. Similar to Fig. 5, except for the Metop rising (a, b) and setting (c, d) occultations, for all latitudes.

Figure 8 .
Figure 8. Similar to Fig. 5, except for the COSMIC L2C setting occultations, for all latitudes.Panels (a, b) and (c, d) correspond to different height intervals and different extrapolation (transition) heights (colors indicated in panels).

Figure A2 .
Figure A2.Atmospheric refractivity (a) and ionospheric electron density (b) used in the modeling.

Figure A3 .
Figure A3.Bending angle profiles obtained from modeling by ray tracing.Black: bending angle obtained with the atmosphere only.Red and green: L1 and L2 bending angles obtained with the atmosphere and ionosphere.Blue: ionospheric-free bending angles.

Figure A4 .
Figure A4.Lower panel: set of ray trajectories between GPS and LEO spanning the troposphere.Upper panel: functions modeling horizontal inhomogeneity of electron density in the ionosphere.

Figure A5 .
Figure A5.The α 1 (h) − α 2 (h) for spherically symmetric (a) and non-spherically symmetric (c1, d1) ionospheric refractivity, for the exponential atmospheric refractivity (black lines) and with added inversion layer (red lines).The errors of the ionospheric correction (c2, d2) in the presence of the inversion layer and horizontally inhomogeneous ionosphere: the standard correction (red lines) and by extrapolation (green lines).