One‐Dimensional Variational Ionospheric Retrieval Using Radio Occultation Bending Angles: 1. Theory

A new one‐dimensional variational (1D‐Var) retrieval method for ionospheric GNSS radio occultation (GNSS‐RO) measurements is described. The forward model implicit in the retrieval calculates the bending angles produced by a one‐dimensional ionospheric electron density profile, modeled with multiple “Vary‐Chap” layers. It is demonstrated that gradient based minimization techniques can be applied to this retrieval problem. The use of ionospheric bending angles is discussed. This approach circumvents the need for Differential Code Bias (DCB) estimates when using the measurements. This new, general retrieval method is applicable to both standard GNSS‐RO retrieval problems, and the truncated geometry of EUMETSAT's Metop Second Generation (Metop‐SG), which will provide GNSS‐RO measurements up to about 600 km above the surface. The climatological a priori information used in the 1D‐Var is effectively a starting point for the 1D‐Var minimization, rather than a strong constraint on the final solution. In this paper the approach has been tested with 143 COSMIC‐1 measurements. We find that the method converges in 135 of the cases, but around 25 of those have high “cost at convergence” values. In the companion paper (Elvidge et al., 2023), a full statistical analysis of the method, using over 10,000 COSMIC‐2 measurements, has been made.


Introduction
The radio occultation instrument on EUMETSAT's Metop Second Generation (Metop-SG) satellites will measure up to 600 km above the Earth's surface, and provide measurements for space weather applications.The retrieval of ionospheric profile information from GNSS radio occultation (GNSS-RO) is well established (e.g., Hajj & Romans, 1998;Schreiner et al., 1999), using the Abel transform, which, for a circularly symmetric electron densities N e and (therefore) refractive indices n, and bending angles α that are available at all impact parameters a, is given by (1) 10.1029/2023SW003572 2 of 15 density profile from topside Incomplete RO data), which iterates toward an electron density profile by modeling the densities above 500 km by a (simplified) VaryChap layer, whose parameters are estimated from Abel-inverted lower altitude electron densities at the previous iteration, is fast enough for operational implementation but prob- ably not yet accurate enough.
In this work, we present a new method based on a variational (or "optimal estimation") approach (Rodgers, 2000).
The one dimensional variational (1D-Var) retrieval approach is used extensively in neutral atmosphere GNSS-RO applications (e.g., Healy & Eyre, 2000;Palmer et al., 2000).Variational retrieval is more flexible than the Abel transform inversion, as it does not rely on an idealized measurement geometry.Specifically, the truncated Metop-SG geometry can be accounted for in a straightforward manner as part of the forward problem,   = () , which maps the electron density profile information, x, to an observation y-bending angles here.This is in contrast to most GNSS-RO ionospheric retrievals, which, with the exception of Hajj and Romans (1998), are based on slant total electron contents (STECs), S.This is defined as the electron density N e integrated along the path P between the GNSS and LEO satellites, thus: However, STEC values are relative quantities.The ionospheric delay experienced by a radio occultation signal depends on its frequency.This would appear to allow the ionospheric effect to be isolated, and the STEC to be calculated, by taking the difference between the delays at two separate frequencies.This is not directly possible, however, because of a remaining a phase ambiguity, which arises because the phase measurements are known only modulo the wavelength (e.g., Dyrud et al. (2008)), and because of differential code biases (DCBs), introduced by both satellite and receiver (e.g., Equation 2 in Montenbruck et al. (2014)).Both these effects need to be accounted for before the STEC can be estimated.In this paper we show how bending angles are related to the vertical derivative of STEC, ∂S/∂a, where a is the impact parameter of the ray path.Taking the derivative removes the sensitivity to any constant offsets (the DCBs) and simplifies the problem.
The standard Abel transform, Equation 1, relates profiles of refractivity and bending angle that range from the surface to infinity.Real observations do not of course extend that far, and bending angles above the maximum height of the Metop-SG measurements, 600 km, are not small enough to be neglected.There are two standard Abel transform methods of generating electron density profiles from such truncated radio occultation signals (e.g., Schreiner et al., 1999).Neither requires absolute TEC, but each faces a difficulty.The first uses bending angles, and requires knowledge of the refractive index at the LEO.The second uses the vertical derivative of the slantwise TEC (or excess phase) (e.g., Equation 3 of Lei et al. (2007)).Unfortunately this necessarily diverges as impact heights approach the LEO, and the (integrable) singularity must be handled somehow.By contrast, in the method proposed in this paper the electron density at any height, including that of the LEO, can be derived from the parameters of the presumed VaryChap layer(s).This also means that electron densities can be inferred above the maximum height of the observations, which makes the method ideally suited for use with truncated bending angles.By construction, the method also enforces positive electron densities, thereby avoiding a failing that can be suffered by Abel inversion in the presence of horizontal gradients (see, e.g., Case 4 in Section 6.1).Finally, the variational method automatically generates measures of the quality of any particular retrieval (number of iterations, cost function at convergence, solution error covariance matrix, etc.)-see Section 2. Note, however, that the retrieval method described in this paper, being one-dimensional, cannot handle horizontal gradients because it assumes spherically symmetric refractivity/electron density fields.This drawback is shared by Abel transform methods.
We note that bending angles are assimilated in most global numerical weather prediction (NWP) data assimilation systems without bias correction, and they are considered "anchor measurements" because they constrain the bias corrections applied to other measurements (e.g., Poli et al., 2010).However, in the ionospheric data assimilation literature, radio occultation STECs with bias corrections (DCBs) are still usually assimilated, for example, Angling and Cannon (2004), Yue et al. (2011), andAngling et al. (2021).
The purpose of this paper is to describe the theoretical basis of the new ionospheric 1D-Var retrieval approach.This includes a description of the multiple Vary-Chap layer electron density profile, the relationship between the derivative of the slant TEC and bending angle, and the 1D-Var assimilation system.A companion paper, Elvidge et al. (2023), presents a full statistical validation of the method.

Data Assimilation Preliminaries
In a variational retrieval system the estimate of the state vector is obtained using the observations combined with a priori information.Since the observations, and a priori information, are only known up to certain levels of confidence these are usually described using probability density functions (PDFs).By assuming Gaussian PDFs, and that the observation and background errors are uncorrelated, Bayes' theorem can be used to show that this problem can be framed as the minimization of a cost function J (e.g., Talagrand, 2010): Here, x b is the background state (specifically, the parameters that define the vertical electron density profile), B is the background error covariance matrix, y o is the set of observations (specifically, bending angles as a function of impact parameter), R is the observation error covariance matrix, and   is the forward operator, which generates pseudo observations from a particular state.
The retrieval solution, x a , is the state which minimizes the cost function, and it should be consistent with both the background x b and the observations y o , to within their expected error statistics.The bending angle profile provides useful information about all the ionospheric variables in the state vector x, and the uncertainty in the estimate of x is significantly reduced as a result of making the measurement.In this context, a useful feature of a 1D-Var retrieval is that it provides the following estimate of the theoretical solution error covariance matrix, A in terms of B and R  = where H is the matrix of partial derivatives of the forward model with respect to the state vector elements:   = ∕ .In well-posed problems the diagonal elements of the A matrix are significantly smaller than the corresponding diagonal elements of the B matrix (A i,i ≪ B i,i ).
The 1D-Var approach should be more robust to measurement noise than the Abel transform because it is a weighted least-squares approach; it will not over-fit the measurements if the assumed observation error statistics, R, are a reasonable approximation of the actual observation error statistics.The retrieval also produces useful diagnostics of the quality of the solution.These include the number of iterations required to converge, and the cost at convergence, J(x a ), the expectation value of which is half the size of the observation vector,  [2 ()] =  ± √ 2 (e.g., Rodgers (2000)), if the assumed error statistics are well specified.
Practical details of the solution method are discussed in Section 5.

Theory
To overcome the need for estimation of the Differential Code Biases (DCBs) in this work the derivative of slantwise TEC S, with respect to the impact parameter a, ∂S/∂a, is assimilated.The latter is the quantity used in the Abel transform solution for refractivity (see Equation 14in Schreiner et al. (1999)), and we will show that, to a good approximation, it is proportional to the difference between the L2 and L1 bending angles, plus a term that involves the electron density at the LEO.
The STEC is the integrated electron density (N e ) along the path P between the GNSS and the LEO satellites (see Equation 2).To a first approximation the delay ϕ i (in m) in the phase of the carrier wave (of frequency f i ), relative to the vacuum, accumulated along the path P is given by where n i is the refractive index at frequency f i , which is approximately given by (e.g., Schreiner et al., 1999) in which the proportionality constant κ ≈ 40.3 m 3 s −2 .
If the electron density is assumed to be only a function of height (i.e., is spherically symmetric), the STEC S between a GNSS satellite at radius r G and a LEO satellite at radius r L is easily shown, from Equation 2 to be given approximately by in which N e (r) is the electron density at radius r, and the impact parameter a is nearly equal to the radial distance to the tangent point, r T .(For example, even for a high electron density of 3 × 10 12 m −3 at a height of 300 km, which is appropriate for the F2 peak in daytime, solar maximum conditions, n i − 1 ≈ − 8 × 10 −5 at the L2 frequency of about 1.2 GHz.This means that a and r T differ by less than 600 m.) Assuming the electron density is zero at the GNSS, integrating Equation 7 by parts gives The derivative of Equation 8with respect to a is: Equation 9 is singular at a = r L , as noted in Lei et al. (2007).This singularity is necessary for the Abel inversion using slantwise TECs to be well-behaved.
Note that the integrals in Equation 9 are closely proportional to the standard Abel transform approximation to the ionospheric bending angle at impact parameter a and frequency f i , α i (a), namely (Angling et al., 2018;Kursinski et al., 1997;Vorob'ev & Krasil'nikova, 1994): 6, n i ≈ 1 (and therefore r T ≈ a) where appropriate, and replacing the infinite upper limits of the integrals by distances to the satellites.
The vertical gradient of the phase delay in Equation 5is ∕ .Therefore, given measurements at two frequencies, f 1 and f 2 , and assuming straight line paths with the same impact parameter, a, we can use Equation 9to write the difference in the phase delay gradients as: There is a fortuitous cancellation of errors relevant to this problem.In the processing of real measurements, bending angle values are derived from the Doppler shift values assuming that the refractive index at the LEO satellite is unity: n(r L ) = 1.For circular orbits this assumption does not affect the impact parameter, a, but it introduces a frequency dependent negative bias, b i , in the observed bending angles.(The proportionality between Doppler shift and impact parameter for circular orbits is also exploited in the Full Spectral Inversion (FSI) technique (Jensen et al., 2003).)This is equal to (see Equation A6 of Schreiner et al. (1999)): which, being inversely proportional to the square of the frequency, cancels out in the usual dual frequency ionospheric correction of the bending angles (Vorob'ev & Krasil'nikova, 1994).The observed ionospheric bending angles,   () , will therefore be related to the true bending angles, α i (a), (see Equation 10), according to Comparison with Equation 11 therefore reveals that Therefore, the ionospheric retrieval can use either observed bending angle differences or the derivative of phase delay differences with respect to impact parameter.

Observations
Accurate estimation of the uncertainty of the observations is crucial to ensure the observations are not over-fitted.
In general it is difficult to estimate the uncertainty of observations, since in most cases there is no additional, independent, reference truth.
As well as errors in the observations themselves, there are also errors arising from the use of a 1D-retrieval rather than a more realistic 3D-retrieval.Such errors obviously depend on the electron density distribution, which varies with location, time of day and solar conditions.To estimate these errors, electron density distributions appropriate to 143 Metop-SG occultations have been generated with 1D and 3D versions of the climatological ionosphere electron density model NeQuick (Nava et al., 2008), at four different solar activity levels, represented by the solar radio flux at 10.7 cm (F 10.7 ).Specifically we use F 10.7 values of 80, 130, 180, 230 solar flux units (sfu) (where 1 sfu = 10 −22 Wm −2 Hz −1 ).
Define   3  (ℎ) as the electron density at height h that results from calculating the phase delays incurred by signals propagating through the ith occultation plane, using 3D NeQuick electron densities, and then passing the delays through an Abel transform to infer an electron density profile at the tangent point.

𝐴𝐴 𝐴𝐴 1𝐷𝐷
(ℎ) is the same, except that it uses a spherically symmetric electron density field, which matches the 1D NeQuick field at the occul tation tangent point.The difference between the two is a measure of the errors in an inferred electron density distribution incurred by the use of a 1D rather than a 3D retrieval.For each altitude h the absolute mean electron density error ϵ across all 4 × 143 = 572 simulated occultations is calculated thus: To ensure that "rare" geometries, such as occultations that cross the day/night terminator, do not distort the statistics, any outliers at each altitude are removed.Outliers are defined as any value which is less than 1.5 times the interquartile range (IQR) below the first quartile (Q1) of the data, or greater than 1.5 times the IQR above the third quartile (Q3), that is, only data, d, in the range are kept at each altitude.The bending angle profile resulting from passing the 1D electron density difference ϵ(h) in Equation 15 through the Abel transform (which is linear in N e ) is shown in Figure 1.The blue vertical line at 2.10 μrad shows the average mean absolute error, and shows that 2.0 μrad is an excellent single value error to use.It is comparable to the value used in neutral atmosphere applications in the middle/upper stratosphere.For example, ECMWF uses 3.0 μrad above about 32 km when assimilating GNSS-RO data.It is clear, however, that the error is height-dependent.Rather than directly using the average errors from Figure 1, which are subject to sampling errors themselves, the following Gaussian fit to the data could be used in a 1D-Var retrieval: bending angle error = 3.8 rad × exp where h is the altitude in km.A minimum value of 1.0 μrad is imposed.Equation 17 is plotted in red in Figure 1.This function provides a similar residual sum of squares error as using a fifth order polynomial fit, but without the associated numerical instabilities.

Electron Density Model
To reconstruct the ionosphere from the observations described in Section 3.2 a model electron density state is needed for the retrieval.Here we assume that the ionosphere can be modeled as a collection of one-dimensional 'Vary-Chap' electron layers.
The standard description (e.g., Wang et al., 2021) of a Vary-Chap electron density profile (described in Reinisch et al. (2007) as a generalization of an α-Chapman profile based on the work by Rishbeth and Garriott (1969)) with a height-dependent scale height, H(h), is given by where In these equations N m is the peak electron density, which occurs at h = h m .H m is the scale height at h m .The estimation of the Vary-Chap parameters in AVHIRO is based on minimizing a least squares cost function, with two terms (Lyu et al., 2019).The first term is the fit to a previous state estimate (rather than a "background" constraint); the second is the fit to the observed STEC measurements, which can be calculated from the L1 minus L2 phase delay measurements.
In practice the height-dependent scale height, H(h), in Equation 18 can be difficult to determine (Kutiev et al., 2009;Nsumei et al., 2012;Wang et al., 2021).Throughout this work the scale height used in the Vary-Chap profile varies linearly with height (above the peak).This means that Equations 18 and 19 reduce to ) , where = log(∕)∕ and ( 21) These equations if h > h m .If h ≤ h m , or if k ≤ 10 −3 , the 'standard' Chapman layer approximation, which is given by the limit as k tends to 0 of the above, applies.This is given by: ) , where is plotted between 100 and 600 km in Figure 2.This one layer Vary-Chap profile provides a good approximation to the standard electron density profile.However, a better approximation can be formed through the addition of multiple Vary-Chap profiles, for example, one for each ionospheric layer (D, E, F1, and F2).It can also be beneficial to introduce a "topside" layer, to try to account for the systematic underestimation of electron densities above the F2 peak height given by the Vary-Chap model.(Prol et al. (2019) show that an F2 scale height that varies quadratically, rather than linearly, with height gives a better fit to topside sounder measurements).Various multi-layer Vary-Chap profiles using the parameters in Table 1 are shown in Figure 3.
The resulting five-layer Vary-Chap model provides an excellent approximation to the ionosphere, with a realistic looking E and F1 region.Using a fifth layer for the topside addresses the underestimation as highlighted by Prol et al. (2019).In practice the addition of the D-region layer has very little impact on the overall profile.Figure 4 shows the difference between a five-layer (F2+F1+E+D+Topside) and four-layer (F2+F1+E+Topside) version, that is, the difference in electron density by including the D-layer.As expected the impact of including a D-Region in the model increases with decreasing altitude, although the maximum difference is approximately 2 × 10 7 m −3 , which is over three orders of magnitude smaller than the absolute density at that altitude.The differences are insignificant, yet they add considerable time to the calculations.(See Section 6.2).

Practical Considerations
The equations of Sections 3.1 and 4 are encoded numerically as follows.The forward model calculates the electron density profile defined by the four Vary-Chap parameters {N m , h m , H m , k} according to Equations 20-25.The integrals needed to calculate the slantwise TEC, S, according to Equation 8 are estimated numerically by Simpson's rule, using at least 30 points between the tangent point and the LEO (counted twice), and 10 points between the LEO and the GNSS (counted once).The integrable singularity at r = a is handled by working in terms of cosh −1 (r/a) rather than r.The procedure is repeated for each Vary-Chap layer.The resulting ∂S/∂a is related to the vertical gradient of the excess phases according to Equation 11, and this in turn is equated to the bending angle differences according to Equation 14 Table 1.
capture the key features of the ionosphere, without complicating the solution by trying to model smaller scale structures at lower altitudes.
The tangent linear model, H, which is needed to minimize the cost function J(x) of Equation 3, and the solution error covariance matrix A via Equation 4, is calculated by evaluating manually differentiated counterparts of the non-linear model.
The cost function J(x) is minimized by means of a standard version of the Levenberg-Marquardt minimization algorithm (Press et al., 1992), which the 'diagonal weighting factor' (usually called λ), which determines how closely the change in x follows the path of steepest descent, is multiplied by 0.1 if the cost function is decreasing, and by 100 if the cost function starts increasing.Iterations are deemed to have converged when changes to the state vector or to the cost function are sufficiently small.In addition, 'unphysical' state vector components, such as negative peak electron densities, are handled during the iterative phase of the minimization, usually by resetting them to a few percent of the associated errors, that is, the square roots of the diagonal elements of B. These are deliberately chosen to be rather large-typically around {5.0 × 10 11 m −3 , 100 km, 20 km, 0.05} for {N m , h m , H m , k} respectively.
The results of this paper have been produced by code written in Python 3.0.Another, widely available version, which is written in Fortran 95, has been part of the Radio Occultation Processing Package (ROPP) since version 11.0 (released January 2022).ROPP is a collection of software, build scripts, test scripts and documentation, which is intended to help users to process RO data.It is maintained, developed and supported by EUMETSAT, through the Radio Occultation Meteorology Satellite Applications Facility (ROM SAF), and freely available to download from its website (ROM SAF, 2023).Users should bear in mind that ROPP is under constant development, and that its code may not therefore exactly match that described in this paper.The differences, however, should be small.

Example Retrievals
The 1D-Var retrieval results are compared with the AVHIRO retrieval (Lyu et al., 2019), and an Abel transform solution.The Abel transform retrievals are not absolute electron density values, because we do not add the electron density at the LEO to the profile.In addition, we have not corrected (calibrated) the bending angles to only include the section of the path below the LEO satellite, by subtracting positive elevation values with the same impact parameter (see Section 3.1 of Schreiner et al. (1999)).However, the Abel solution should indicate whether or not the 1D-Var results look reasonable.Our implementation of the Abel transform assumes that the L2-L1 bending angle differences vary linearly in the vertical between observations.This means that the Abel transform is linear, so the retrieved electron density profile can be computed as a matrix multiplied by a vector of (differenced) bending angles.
The performance of the 1D-Var is illustrated with four cases (see Elvidge et al. (2023) for the complete analysis), which have been chosen for their differing retrieval characteristics and convergence properties.

Case 1
The first case, shown in Figure 5, can be considered a good retrieval, and is the same example as shown in the top right panel of Figure 4 in Lyu et al. (2019).The 1D-Var retrieval uses two layers, an F2 layer and F1 layer, two layers being a reasonable compromise between accuracy (which increases with number of layers) and robustness and CPU time (which decrease with it)-see Section 6.2.(Also see Elvidge et al. (2023), which contains an extensive study of the impact of adding more and more layers to the 1D-Var retrieval method described in this paper).A one-off Abel Transform solution has also been added, to replicate Figure 4 of Lyu et al. (2019) more completely.The three approaches agree very similarly around the profile peak (between approximately 200 and 400 km) but the 1D-Var solution differs from the AVHIRO and the Abel transform solutions outside this region.

Case 2
The second case, shown in Figure 6, compares the 1D-Var retrieval with those from COSMIC-2 (UCAR Cosmic Program, 2019) and a closely located ionosonde.Ionosondes step through a range of HF frequencies transmitted vertically upwards and measure the return echoes, which enables them to image the vertical profile of the ionosphere up to the peak density.(Ionosonde density profiles above the peak are less reliable, because they depend on a background model, which is why they are plotted as dashed lines in this paper.)Ionosonde observations are widely used as reference observations for comparative studies (e.g., Feltens et al., 2011;Elvidge et al.. 2014;Scherliess et al., 2011).Here an ionosonde observation (profile) from within 200 km of the location of the occultation have been used to demonstrate the 1D-Var retrieval (A full statistical study is undertaken in Elvidge et al. ( 2023)).Also shown is the convergence in bending angle space, that is, the observed (blue) and forward modeled solution (green) bending angle differences α 2 − α 1 .
In this case the 1D Var solution is in very close agreement with the ionosonde observations throughout the whole profile.Above approximately 300 km and below 160 km there are some deviations from the observations, but these are small.The COSMIC-2 retrieval is also very good, but overestimates the peak electron density.The COSMIC-2 profile shows more structure below 150 km, presumably resulting from the variability in the α 2 − α 1 observations in that region.The 1D-Var solution is the result of fitting observations between 500 and 175 km (see Section 5), and therefore has no knowledge of this structure.Even so, the forward modeled bending angle differences give a reasonable fit to the observations.

Case 3
The third case, Figure 7, again compares the 1D-Var retrieval with COSMIC-2 and a nearby ionosonde (in the same manner as described in Section 6.1.2) and both the COSMIC-2 and 1D-Var retrievals are fairly similar.The peak heights of the profiles, at about 200 km, are similar and both are slightly higher than that reported by the ionosonde.(But note that autoscaled ionosonde observations should be treated carefully, and the errors in the F2 peak height can be reasonably large (as much as 10%, according to Themens et al. (2022)).The COSMIC-2 retrieval slightly overestimates the peak density whilst the 1D-Var slightly underestimates it.Between the peak of the electron density profile and about 150 km, the COSMIC-2 retrieval more closely matches the observations.Above about 275 km the two retrievals are very similar.(Differences from the ionosonde profile in this region should be treated skeptically, for reasons explained in Section 6.1.2).

Case 4
In the final case, shown in Figure 8, the COSMIC-2 retrieval is generally closer to the ionosonde observation than the 1D-Var retrieval.Again, both retrievals have similar heights of peak density, but the COSMIC profile more closely matches the observations above and below the peak.However the COSMIC-2 retrieval returns negative electron densities below 100 km, whereas the 1D-Var retrieval cannot, for reasons described in Section 5.

Convergence and Computational Expense
Elvidge et al. ( 2023) provides a full statistical analysis of the results of the 1D-Var retrievals compared to Abel-transform-based COSMIC-2 retrievals.This section simply provides a flavor of the 1D-Var ionospheric retrieval method's performance on 143 bending angle profiles provided by the Institut D'Estudis Espacials de Catalunya (IEEC) in Barcelona, Spain.These are COSMIC-1 measurements from 18 September 2011.The data contains the "geometry free" (ϕ 1 − ϕ 2 ) phase differences as function of impact parameter, a, up to the COSMIC-1 satellite, which operated at an altitude of about 800 km.
The simulated differenced bending angles α 2 − α 1 are computed with a finite difference approximation to ∂(ϕ 1 − ϕ 2 )/∂a.The vertical separation between the bending angles is typically 500 m.A common radius of curvature of 6,371 km is used for all cases since it was not included in the test data sets.
Although adding layers seems to improve the overall shape of the retrieved ionosphere, this comes at the cost of higher computational time and number of iterations required for the 1D-Var to converge.This is shown in Table 2.Note that the maximum number of iterations was set to 50.Any retrievals needing more iterations than this were recorded as convergence failures, and discarded from further analysis.On the basis of the 143 profiles studied here, two layers generate a reasonable fit to the observations at a moderate computational cost, which is suitable for illustrating the method.In the companion paper Elvidge et al. (2023), however, which uses a much better background model on a much larger sample, four layers (F2, F1, E and topside) work better, having a smaller RMS error with respect to ionosondes (albeit at a lower convergence rate).The correct number to use in any particular application will depend on the balance of the competing demands of speed, robustness and accuracy that the user feels appropriate for that application.
It can be seen from Table 2 that the average time per iteration increases as the number of layers increases.
However after an immediate jump in the average number of iterations when using one-and two-layer models, the average number of iterations for converging observations remains fairly steady (at around 30).This leads to a monotonic increase in the average total length of time needed for an observation to converge using the 1D-Var, as shown in Figure 9.In addition, the percentage of observations that converge continually decreases as the number of layers increases, from 98.6% with just one layer to only 58.7% with five layers.More realistic retrieved electron density profiles therefore come at a cost in CPU computational time and robustness.

Discussion and Conclusions
This paper has described a new 1D-Var ionospheric retrieval method.This approach can be applied to the Metop-SG measurement geometry, which will truncate the ionospheric GNSS-RO measurements around 600 km above the surface, as well as 'standard' GNSS-RO measurement geometry, in which measurements are assumed to be available at all heights.
The original plan was to use the difference between the slantwise TECs of the L1 and L2 radio signals, expressed as a function of impact parameter, but the L2-L1 observed bending angle differences have been used instead.This approach is closer to the way that GNSS-RO data is used in neutral atmosphere applications, such as operational NWP.
A new 1D ionospheric bending angle forward operator has been developed, which computes L2-L1 bending angle differences as a function of impact parameter.This operator assumes that the ionospheric electron density can be modeled as a collection of 'Vary-Chap' layers.As suggested by a comparison of bending angles generated by 1D-and 3D-NeQuick electron density models, the L2-L1 bending angle uncertainty is assumed to be a constant 2.0 μrad, and vertical error correlations are neglected.Some experiments were made using the more complicated error formula of Equation 17, but the results were found to differ little from the same experiments using 2.0 μrad.For simplicity and reproducibility, therefore, and to illustrate the method, the constant value was adopted.
We find that gradient-based minimization techniques can be successfully applied to this retrieval problem, and that the non-linear, nested exponential nature of the Vary-Chap profiles (see Equation 20) does not cause problems.This is in contrast to the assertion of Lyu et al. (2019).Typically, at least two out of three retrievals converge within 50 iterations, although this convergence rate decreases as more Vary-Chap layers are introduced.The iteration count, and the high final cost functions of the retrievals that fail to converge, provide useful diagnostic information, as does the automatically generated solution error covariance matrix A.
A few example 2-layer 1D-Var retrievals have been compared against the results of the AVHIRO method, Abel transform retrievals, and nearby ionosondes.In general the 1D-Var method performs at least as well as AVHIRO and the Abel transform, and produces a close fit to the observed bending angle differences α 2 − α 1 in the region where it is supposed to.By construction, 1D-Var also avoids one drawback of the Abel inversion technique, namely the generation of negative electron densities.Adding more layers can improve the fit, but at the cost of longer CPU times and a lower convergence rate.
The a priori information is essentially a first guess used to start the minimization, rather than a strong constraint on the final 1D-Var solution, because the background errors are (deliberately) rather large.Better a priori information might speed up convergence but would not necessarily improve the 1D-Var solution.It would, however, make it easier to screen out poor bending angles at the start of the minimization.
The key suggestion of this paper is that it may be possible to use differenced (L1-L2) bending angles in ionospheric data assimilation (DA) systems.Usually, such systems assimilate slantwise TEC values, but these require a correction for the Differential Code Biases (DCBs).Constant DCBs during the occultation, however, will not affect the bending angles derived from raw phase delays, and can therefore be ignored.This suggestion should at some point be tested in the context of a more formal data assimilation system.

Figure 1 .
Figure 1.Mean absolute difference in bending angles derived from electron density differences ϵ(h) defined by Equation 15.Vertical average in blue.Fitting function in Equation 17 in red.

Figure 2 .
Figure 2. of a Vary-Chap profile with "standard" parameters given by Equation 26.

Figure 3 .
Figure 3. Examples of multi-layer Vary-Chap profiles with the parameters given inTable 1.

Figure 6 .
Figure 6.Comparison between the 1D-Var electron density profile retrieval with the Abel transform retrieval of a COSMIC-2 profile and a nearby ionosonde.Also shown are the observations and the forward modeled solution in bending angle space (i.e., α 2 − α 1 ).(Ionosonde profiles above the peak appear as dashed lines, as they are less reliable there).

Figure 8 .
Figure 8.Comparison between the 1D-Var electron density profile retrieval with the Abel transform retrieval of a COSMIC-2 profile and a nearby ionosonde.(Ionosonde profiles above the peak appear as dashed lines, as they are less reliable there).

Figure 7 .
Figure 7.Comparison between the 1D-Var electron density profile retrieval with the Abel transform retrieval of a COSMIC-2 profile and a nearby ionosonde.(Ionosonde profiles above the peak appear as dashed lines, as they are less reliable there).

Table 2
Average Time Per Iteration, Percentage of Observations Which Converge, and Statistics of the Number of Iterations Needed for Convergence, for up to Five Vary-Chap Layers Figure 9. Average time taken per iteration (blue) and average length of time for an observation to converge (red).