Variational Assimilation of Radio Occultation Observations into Numerical Weather Prediction Models: Equations, Strategies, and Algorithms

Abstract

We review different approaches to the variational assimilation of radio occultation (RO) observations into models of global atmospheric circulation. We derive the general equation for the bending angle that reduces to the Abel integral for a spherically layered atmosphere. We review the full 3-D observation operator for bending angles, which provides the strictest solution, but is also most computationally expensive. Commonly used is the 2-D approximation that allows treating rays as plane curve. We discuss a simple 1-D approach to the assimilation of bending angles. The observation operator based on the standard form of the Abel integral has a disadvantage, because it cannot account for waveguides. Alternative approaches use 1-D ray-tracing. The most straightforward way is to use the same framework as for the 3-D observation operator, with the refractivity field reduced to a single profile independent from the horizontal coordinates. An alternative 1-D ray-tracing approach uses the form of ray equation in a spherically layered medium that uses an invariant. The assimilation of refractivity has also 1-D and 3-D options. We derive a new simple form of the refractivity-mapping operator. We present the results of numerical tests of different 3-D and 1-D observation operators, based on Spire data.

1. Introduction

Radio occultation observations have always been looked at as an important source of data for numerical weather prediction (NWP) [1,2,3,4,5]. A general approach to the use of observations in NWP is based on the variational data assimilation [6]. A way of the implementation of this approach for radio occultation (RO) observation was first discussed by Eyre [7], who applied the finite difference technique for the explicit derivation of the 3-D adjoint operator based on the ray-tracing.

Zou et al. [8] described assimilation of RO data into NWP model using the simplest observation operator: the refractivity profiles retrieved in the assumption of the spherically symmetric atmosphere were assumed to be the observables. Zou et al. [9] developed a 3-D observation operator based on the ray-tracing, and its adjoint. The bending angle was defined as the angle between the initial and final ray direction. This operator was successfully tested by Zou et al. [10] and then used by Liu et al. [11] and enhanced by Liu and Zou [12]. This approach is also referred to as 2-D approach, because it is a good approximation to consider rays as plane curves, i.e., to neglect the transverse displacement of a ray.

Kuo et al. [13] discussed possible RO data assimilation strategies. The following candidates for observables were discussed: (1) raw phases and amplitudes in the framework of the diffraction theory; (2) raw phases in the framework of the geometrical optics (GO) under the assumption of single-ray propagation; (3) L1 and L2 bending angles; (4) ionospheric-corrected bending angle; (5) retrieved refractivity; and (6) retrieved temperature. A conclusion was made that the assimilation of bending angles, as was done by Zou et al. [10] should be preferred.

Gorbunov and Kornblueh [14] developed and implemented the strict 3-D observation operator and its adjoint. The operator is also based on ray-tracing. However, the bending angle and impact parameter are defined in terms of Doppler frequency. These are referred to as “effective” bending angle and impact parameter. The reason is that in a generic 3-D atmosphere, the ray impact parameter is not invariant for each ray, and, therefore, the impact parameters at the ray final and end points are different and cannot be derived from the observed Doppler frequency without additional assumptions. We define “effective” quantities as those that would correspond to the observed Doppler frequency, if the atmosphere were spherically symmetric. This does not involve any assumptions, because this is the definition of the observables, which is implemented in the same way in the data processing and in the observation operator and, therefore, does not result in any mismodeling. The code developed by Gorbunov and Kornblueh [14] was later modified and included into the operational data assimilation system of German Weather Service [15,16]. Although the 3-D approach is most accurate, its obvious disadvantage consists in its high computational demands.

Kursinski and Hajj [17], Healy et al. [18] followed the 1-D approach based on the retrieved refractivity. Von Engeln et al. [19] introduced a 1-D approach based on bending angles that are linked to the refractivity profile by the Abel transform. This approach was adopted by Healy and Thépaut [20], who used an approximation of the Abel integral based on piecewise exponential refractivity profile. Cucurull et al. [21,22,23] also used 1-D bending angle operator, but with a different numerical scheme of the evaluation of the Abel integral, based on the coordinate change eliminating the pole of the integral kernel. Burrows et al. [24] developed improvement of the interpolation scheme to reduce numerically induced biases.

Syndergaard et al. [25] introduced a refractive index mapping operator, which maps the 3-D field of model refractivity to the retrieved 1-D vertical refractivity profile. The relation between the 3-D field and 1-D profile involves the forward and inverse parts, the former representing the integration of refractivity along the rays, and the letter being the Abel inversion. Using an approximate form of ray trajectories, it was shown that it is possible to evaluate the matrix of the operator once and then use it without additional computational expenses. Sokolovskiy et al. [26,27] performed numerical simulations to test this approach with some modifications: phase observable was used instead of the retrieved refractivity. The phase observable is defined as the integral of the 1-D Abel-retrieved refractivity along straight lines. Sokolovskiy et al. [26] suggested terms “local” for 1-D operators and “non-local” for 3-D and 2-D operators. A further comparison of local and non-local operators was performed by Ma et al. [28], who confirmed that this approach allows a computationally inexpensive account of the horizontal structure using the non-local phase excess operator. Liu et al. [29] evaluated a non-local quasi-phase operator in assimilation of the CHAMP GPS refractivity using the Weather Research and Forecasting (WRF) ensemble data assimilation system. A conclusion was made that the use of the non-local RO quasi-phase operator can significantly improve the assimilation of RO data in the WRF ensemble assimilation system in the middle and lower troposphere. A further development of the mapping approach was made by Aparicio [30] who introduced the idea of using a linear or quadratic form based on the regression in order to represent the refractivity field in the vicinity of the perigee.

Healy et al. [31] investigated and compared two approaches to the 2-D bending angle observation operator. The first approach is based on the approximate form of ray and bending angle equation represented as the Abel integral with varying impact parameter. The second approach uses the accurate ray equations in the polar coordinates. It was demonstrated that as compared to 1-D approach, 2-D observation operator has smaller modeling errors.

2. Basic Equations

2.1. Ray Equations

The electromagnetic field in a medium is described by the Maxwell equations system. Using simple approximations, it is possible to reduce the system to the following scalar (Helmholtz) equation [32]:

$$\mathcal{H}\left(\mathbf{x},\widehat{\mathbf{p}}\right)u=0,$$

(1)

where u is any component of the electric field, $\mathcal{H}$ is the Hamilton function, $\mathbf{x}=\left(x^i\right)$ is the vector of the Cartesian coordinates, and $\widehat{\mathbf{p}}=\left({\widehat{p}}_i\right)$ is the covector of the momentum operators [33]:

$${\widehat{p}}_i=\frac{1}{\widehat{i}k}\frac{\partial }{\partial x^i},$$

(2)

where $\widehat{i}$ stays for the imaginary unit, $k=2\pi /\lambda$ is the wavenumber, and $\lambda$ is the wavelength. The Hamilton function has the following form:

$$\mathcal{H}\left(\mathbf{x},\mathbf{p}\right)=\frac{1}{2}\left({\mathbf{p}}^2-n^2\left(\mathbf{x}\right)\right),$$

(3)

where $n\left(\mathbf{x}\right)$ is the refractivity field. The solution of the Helmholtz equation has the form of $u\left(\mathbf{x}\right)=A\left(\mathbf{x}\right)\exp \left(ik\mathrm{\Psi }\left(\mathbf{x}\right)\right)$, where $A\left(\mathbf{x}\right)$ is the amplitude, and $\mathrm{\Psi }\left(\mathbf{x}\right)$ is the eikonal. The lowest order asymptotic results in the eikonal equation:

$$\mathcal{H}\left(\mathbf{x},\nabla \mathrm{\Psi }\right)=0,n^2={\left(\nabla \mathrm{\Psi }\right)}^2.$$

(4)

This a partial differential equation, whose solution can be obtained from the characteristic equations, which form the Hamilton system:

$$\begin{array}{cc}\hfill & \dot{\mathbf{x}}=\frac{\partial \mathcal{H}}{\partial \mathbf{p}},\dot{\mathbf{p}}=-\frac{\partial \mathcal{H}}{\partial \mathbf{x}},\dot{\mathrm{\Psi }}=\mathbf{p}\dot{\mathbf{x}},\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill & \dot{\mathbf{x}}=\mathbf{p},\dot{\mathbf{p}}=n\nabla n,\dot{\mathrm{\Psi }}=n^2,\hfill \end{array}$$

(6)

where $\mathbf{p}$ is the classical momentum. Because $\left|\mathbf{p}\right|=\left|\nabla \mathrm{\Psi }\right|=n$, as follows from (4) and (6), we arrive at the following differential relation between the parameter t of system (6), the ray arc length s, and the eikonal:

$$dt=\frac{ds}{n},d\mathrm{\Psi }=nds.$$

(7)

The equation for the amplitude follows from the ray equations and the energy conservation, which, in the geometric optical approximation, requires that the energy should be transferred along rays.

Equation (6) has a form that is specific for the Cartesian coordinates. Consider an arbitrary coordinate system with the metric tensor $g_{ij}$: ${ds}^2={dx}^i{g}_{ij}{dx}^j$, where we, according to the rules of the tensor calculus, follow the Einstein notation implying the summation over each pair of upper and lower indexes of the same name. We define the momentum by the relation $p_i=g_{ij}{\dot{x}}^j$. Because form $\mathbf{p}d\mathbf{x}$ is invariant with respect to a coordinate change, the transform to the new coordinates $\left(\mathbf{p},\mathbf{x}\right)$ is canonical, and the canonical form of the Hamilton system also remains invariant [34], provided that the Hamilton function is defined as follows:

$$\mathcal{H}\left(\mathbf{x},\mathbf{p}\right)=\frac{1}{2}\left(p_i{g}^{ij}{p}_j-n^2\left(\mathbf{x}\right)\right),$$

(8)

where $g^{ij}$ is the matrix inverse to $g_{ij}$. This results in the following form of the ray equations:

$$\begin{array}{cc}\hfill & {\dot{x}}^i=\frac{\partial \mathcal{H}}{\partial p_i}=g^{ij}{p}_j,\hfill \\ \hfill & {\dot{p}}_i=-\frac{\partial \mathcal{H}}{\partial x^i}=n\frac{\partial n}{\partial x^i}-\frac{1}{2}{p}_k\frac{\partial g^{kj}}{\partial x^i}{p}_j.\hfill \end{array}$$

(9)

The 2-D approximation [35] allows treating rays as plane curves. Consider polar coordinates $\left(r,\theta \right)$ with the metric tensor:

$$g_{ij}=\left(\begin{array}{cc}1& 0\\ 0& r^2\end{array}\right),g^{ij}=\left(\begin{array}{cc}1& 0\\ 0& r^{-2}\end{array}\right).$$

(10)

Then we have the following equations:

$$\begin{array}{cc}\hfill & p_{\theta }=r^2\dot{\theta }=nr\frac{rd\theta }{ds}=nr\sin \psi ,\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill & \dot{p_{\theta }}=n\frac{\partial n}{\partial \theta },\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill & {\dot{p}}_r=\ddot{r}=n\frac{\partial n}{\partial r}+\frac{p^2}{r^3}.\hfill \end{array}$$

(13)

where $\psi$ is the angle between ray direction $\dot{\mathbf{x}}$ and the local vertical $\mathbf{x}$. The angular component of the momentum $p_{\theta }$ coincides with the ray impact parameter p, which is invariant in a spherically layered medium. The equation for $p_{\theta }$ does not include ${\dot{p}}_r$ and can be integrated separately. Using Equation (11) and relation $d{s}^2=d{r}^2+r^{2}d{\theta }^2$, we arrive at the following equation:

$$\frac{dr}{d\theta }=\pm \frac{r}{p}\sqrt{n^2{r}^2-p^2}.$$

(14)

Hereinafter, the upper sign relates to the ascending part of the ray trajectory from the perigee to the receiver, and the lower sign relates to descending part of the ray trajectory from the transmitter to the perigee. In the general case of a medium with horizontal gradients, this equation should be complemented with the dynamic equation for p. Equation (12) describes the variations of the ray impact parameter caused by horizontal gradients of refractivity [36,37,38].

In the perigee point, $\psi =\pi /2$ and, therefore, $p=nr$ and from Equation (13) it follows that

$$\ddot{r}=n\frac{\partial n}{\partial r}+\frac{n^2}{r}=\frac{n}{r}\frac{dx}{dr},$$

(15)

where $x=nr$ is the refractive radius.

2.2. Bending Angle

Using the ray equations obtained above, it is straightforward to arrive at the bending angle expression. We will start with the differential equations for the bending angle:

$$\begin{array}{cc}\hfill & dϵ=d\left(\psi +\theta \right),\psi =\arcsin \frac{p}{nr},\hfill \\ \hfill & dϵ=-d\left(\psi +\theta \right),\psi =\pi -\arcsin \frac{p}{nr},\hfill \end{array}$$

(16)

written for the ascending and descending part of the ray trajectory, respectively. These equations can be combined as follows:

$$dϵ=d\arcsin \frac{p}{nr}\pm {\left(\frac{dr}{d\theta }\right)}^{-1}dr.$$

(17)

Taking into account Equation (14), we obtain the following equation:

$$\begin{array}{cc}\hfill dϵ& =\frac{nr}{\sqrt{n^2{r}^2-p^2}}\left(\frac{dp}{nr}-\frac{p}{n{r}^2}dr-\frac{p}{n^{2}r}\left(\frac{\partial n}{\partial r}dr+\frac{\partial n}{\partial \theta }d\theta \right)\right)+\frac{pdr}{r\sqrt{n^2{r}^2-p^2}}\hfill \\ \hfill & =\frac{1}{\sqrt{n^2{r}^2-p^2}}\left(dp-\frac{p}{n}\left(\frac{\partial n}{\partial r}dr+\frac{\partial n}{\partial \theta }d\theta \right)\right).\hfill \end{array}$$

(18)

Together with the dynamic equation for impact parameter (12), this can be transformed as follows:

$$\begin{array}{cc}\hfill dp& =\frac{\partial n}{\partial \theta }\frac{ds}{d\theta }d\theta =\frac{\partial n}{\partial \theta }\sqrt{{\left(\frac{dr}{d\theta }\right)}^2+r^2}d\theta =\frac{\partial n}{\partial \theta }\frac{n{r}^2}{p}d\theta ,\hfill \\ \hfill dϵ& =\frac{1}{\sqrt{n^2{r}^2-p^2}}\left(\frac{\partial n}{\partial \theta }\frac{n{r}^2}{p}d\theta -\frac{p}{n}\left(\frac{\partial n}{\partial r}dr+\frac{\partial n}{\partial \theta }d\theta \right)\right).\hfill \end{array}$$

(19)

In addition, taking into account Equation (14) once again, we arrive at the final expression:

$$dϵ=\left(-\frac{p}{n}\frac{\partial n}{\partial r}\pm \frac{\partial n}{r\partial \theta }\frac{\sqrt{n^2{r}^2-p^2}}{n}\right)\frac{dr}{\sqrt{n^2{r}^2-p^2}}.$$

(20)

At an arbitrary trajectory point, we introduce local Cartesian coordinates $\left(dr,rd\theta \right)$. The ray direction vector in this basis equals $\left(\cos \psi ,\sin \psi \right)$. The expression in the brackets in Equation (20) equals the scalar product of covector $\nabla n$ and vector

$$\left(-\frac{p}{n},\pm \frac{\sqrt{n^2{r}^2-p^2}}{n}\right)=r\left(-\frac{p}{rn},\pm \sqrt{1-{\left(\frac{p}{rn}\right)}^2}\right)=r\left(-\sin \psi ,\cos \psi \right),$$

(21)

which is normal to the ray. In a spherically layered medium, this expression reduces to the standard formulas for the bending angle in terms of the integration over r and over x:

$$\begin{array}{cc}\hfill ϵ\left(p\right)& =-2p\underset{r_0}{\overset{\infty }{\int }}\frac{d\ln n}{dr}\frac{dr}{\sqrt{n^2{r}^2-p^2}}.\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill ϵ\left(p\right)& =-2p\underset{p}{\overset{\infty }{\int }}\frac{d\ln n}{dx}\frac{dx}{\sqrt{x^2-p^2}}.\hfill \end{array}$$

(23)

2.3. Influence of Waveguides

The dependence of the refractive radius $x\left(r\right)=n\left(r\right)r$ from the geometrical radius r provides a comprehensive description of refraction in a spherically layered medium. As indicated by Equation (15), the derivative $dx/dr$ at the ray perigee point, where $\dot{r}=0$, allows distinguishing three situations: (1) $\ddot{r}\sim dx/dr>0$ – normal refraction (ray ascends); (2) $\ddot{r}\sim dx/dr=0$ – critical refraction (ray slides along the spherical surface); (3) $\ddot{r}\sim dx/dr=0$ – super-refraction (ray descends).

Because asymptotically for large height, $n\left(r\right)$ decreases exponentially, and the refractive radius becomes close to the geometrical radius, we can point out that the situation with waveguides and super-refraction corresponds to non-monotonic profiles of $x\left(r\right)$. For each ray, by virtue of Equation (11), $x\ge p$. This leads to the consideration of super-refraction illustrated in Figure 1. In some situations, to a single value of the impact parameter, multiple rays may correspond. However, only one of them can be observed from space, while the other rays are trapped by the waveguide.

Note, the above derivation of the bending angle expression (22) only relied upon the fact that the ray has a single perigee point and reaches the upper border of the atmosphere. Anyway, rays with multiple perigee points must be trapped inside a waveguide and have an infinite bending angle. This expression can, therefore, be applied for the ray with impact parameter $p_2$ in Figure 1. The application of expression (23) in this situation is, however, not straightforward, because x is not a unique coordinate and $n\left(x\right)$ is a multi-valued function. This integral can still be understood as a first-kind integral over a manifold, where $dx$ can change its sign. In this sense, the bending angle expression can be rewritten as follows:

$$ϵ\left(p\right)=-2p\underset{p}{\overset{\infty }{\int }}\frac{d\ln n}{\sqrt{n^2{r}^2-p^2}}.$$

(24)

In solving the inverse problem, waveguide must lead to a systematic negative refractivity error [39]. To see this, consider a waveguide in the radius interval from $r_1$ to $r_2$. As follows from the above discussion, $n\left(r_1\right)r_1=n\left(r_2\right)r_2=x_s$. Under realistic conditions, refractivity profiles monotonically decrease with height. Assuming that $r_1<r_2$, we see that $n\left(r_1\right)>n\left(r_2\right)$.

Consider profile $n^{*}\left(x\right)$, composed of profiles $n\left(x\right)$ below and above the waveguide. At the point $x_s$, the logarithm of this profile has a negative jump $\Delta \ln n=\ln n\left(r_2\right)-\ln n\left(r_1\right)$. To accurately reconstruct profile $n^{*}\left(x\right)$, the following bending angle profile is needed:

$${ϵ}^{*}\left(p\right)=\left\{\begin{array}{cc}\hfill & ϵ\left(p\right),p>x_s\hfill \\ \hfill & -2p\underset{p}{\overset{x_s}{\int }}+\underset{x_s}{\overset{\infty }{\int }}\frac{d\ln n^{*}}{dx}\frac{dx}{\sqrt{x^2-p^2}}+\frac{2p\left|\Delta \ln n\right|}{\sqrt{x_s^2-p^2}},p<x_s\hfill \end{array}\right.$$

(25)

where the additional term accounts the delta-function in the derivative of $\ln n$ for the rays with perigee points below the waveguide. The actual observed bending angle profile for $p<x_s$ can be represented as follows:

$$ϵ\left(p\right)=-2p\underset{p}{\overset{x_s}{\int }}+\underset{x_s}{\overset{\infty }{\int }}\frac{d\ln n^{*}}{dx}\frac{dx}{\sqrt{x^2-p^2}}+2p\underset{n_1}{\overset{n_2}{\int }}\frac{\left|d\ln n\right|}{\sqrt{x^2-p^2}},$$

(26)

where the last term relates to the piece of the profile n inside the waveguide, where $x>x_s$. Therefore, we arrive at the following estimate:

$$2p\underset{n_1}{\overset{n_2}{\int }}\frac{\left|d\ln n\right|}{\sqrt{x^2-p^2}}<\frac{2p\left|\Delta \ln n\right|}{\sqrt{x_s^2-p^2}}.$$

(27)

This indicates that $ϵ\left(p\right)<{ϵ}^{*}\left(p\right)$ for $p<x_s$. This leads to the fact that in the presence of waveguide the Abel inversion of bending angle profiles will result in negative retrieval error of $n\left(r\right)$, while the waveguide structure cannot be retrieved.

3. Observation Operators

3.1. 3-D Assimilation of Bending Angles

3.1.1. The Physical Model of RO Observations

The physical model of RO observations is based on the concept of electromagnetic waves propagating through the inhomogeneous atmosphere. However, it is not expedient to take the observed amplitudes and excess phases as observable and construct the observation operator based on the wave equation in an inhomogeneous medium. Such a formulation would be too detailed, too non-linear, and too unstable with respect to small perturbations. A much better strategy is based on the concept of the ray structure of the wave field that can be very efficiently retrieved using methods based on Fourier Integral Operators [40,41,42,43,44,45,46]. All these methods are closely related to each other and represent different approximations or different coordinate representations of the same solution [47]. The highest possible vertical resolution of these methods is even higher than that required by numerical weather prediction systems. Currently, the retrieval of neutral atmospheric bending angle profile is a commonly accepted procedure of the preparation of RO observation for its further assimilation. Its first step is the retrieval of bending angle profiles in the two channels, followed by the ionospheric correction [48], statistical regularization and noise reduction [49,50,51,52,53,54,55].

A common observation is that, due to a complicated architecture of RO experiments, any formulation of observation operator involves approximations. Given retrieved neutral bending angle, it is sufficient to formulate the forward model in terms of the relative Doppler frequency shift expressed as a function of satellite orbit data and ray geometry:

$$d=\frac{{\omega }_R-{\omega }_T}{{\omega }_T}=\frac{{\mathbf{V}}_{\mathrm{T}}·{\mathbf{u}}_{\mathrm{T}}-{\mathbf{V}}_{\mathrm{R}}·{\mathbf{u}}_{\mathrm{R}}}{c},$$

(28)

where ${\omega }_R$ is the frequency of the received signal, and ${\omega }_T$ is the frequency of the transmitted signal, ${\mathbf{V}}_{\mathrm{T},\mathrm{R}}$ are the transmitter and receiver velocities, ${\mathbf{u}}_{\mathrm{T},\mathrm{R}}$ are unit vectors of the ray directions at the transmitter and receiver, c is the light velocity in a vacuum.

The first numerical simulations of RO sounding of the Earth’s atmosphere [1] involved the iterative solution of the 3-D ray boundary problem, referred to as the shooting method of locating the ray connecting two prescribed satellite positions. This led to the original formulation of the observation operator based on the excess phase $\mathrm{\Psi }\left(t\right)$ or its normalized derivative $d\left(t\right)=(1/c)d\mathrm{\Psi }\left(t\right)/dt$ as functions of time, together with the satellite orbit data ${\mathbf{x}}_{\mathrm{T},\mathrm{R}}\left(t\right)$. This was thought to allow the most precise geometrical parameterization of the assimilation problem. However, the analysis of the GPS/MET RO data [56,57] indicated that multipath propagation effects play a significant role in the troposphere. This aggravates the solution of the boundary problem: (1) the model must be capable of locating multiple rays, whose number is unknown a priori; (2) there is the structural instability: at some points even small perturbations of the atmospheric state would result in the change of the number of rays; (3) the measured Doppler shift or excess phase cannot be used directly because the phase of the wave field in multipath zones is a strongly non-linear combinations of phases of the multiple interfering rays, which need to be separated (cf. the above discussion of the observation operator based on the wave equation); (4) the measured Doppler shift and phase excess need the ionospheric correction, while the most efficient ionospheric correction is the linear combination of bending angles with the same impact parameter [48]. Finally, the solution of the multipath problem [40,41,42,43,44,45,46] indicated that the optimal observable to be assimilated is the bending angle as a function of the impact parameter, because the latter allows both unfolding the multipath propagation effects and reducing the ionospheric correction errors.

The bending angle is the angle between ${\mathbf{u}}_{\mathrm{T}}$ and ${\mathbf{u}}_{\mathrm{R}}$, and, in a general case, it cannot be determined from Equation (28), because two unit vectors have four degrees of freedom, while we have just one scalar equation. However, in the observation operator we can mimic the standard processing scheme based on the sphericity assumption [14]. In a spherically layered atmosphere, ray is a plane curve and impact parameters at the transmitter and receiver are equal to each other. Therefore, we introduce the “effective” ray directions ${\tilde{\mathbf{u}}}_{\mathrm{T},\mathrm{R}}$ restricted to be unit vectors and satisfying Equation (28) and the above sphericity assumption:

$$\begin{array}{cc}\hfill & \frac{{\mathbf{V}}_{\mathrm{T}}·{\tilde{\mathbf{u}}}_{\mathrm{T}}-{\mathbf{V}}_{\mathrm{R}}·{\tilde{\mathbf{u}}}_{\mathrm{R}}}{c}=d,\hfill \\ \hfill & {\tilde{\mathbf{u}}}_{\mathrm{T}}·{\tilde{\mathbf{u}}}_{\mathrm{T}}={\tilde{\mathbf{u}}}_{\mathrm{R}}·{\tilde{\mathbf{u}}}_{\mathrm{R}}=1,\hfill \\ \hfill & {\mathbf{x}}_{\mathrm{R}}\times {\tilde{\mathbf{u}}}_{\mathrm{R}}={\mathbf{x}}_{\mathrm{T}}\times {\tilde{\mathbf{u}}}_{\mathrm{T}},\hfill \end{array}$$

(29)

where ${\mathbf{x}}_{\mathrm{T},\mathrm{R}}$ are the satellite coordinates with respect to the local curvature center of the Earth’s surface cross-section by the occultation plane [58,59]. This system provides six equations defining six components of vectors ${\tilde{\mathbf{u}}}_{\mathrm{T},\mathrm{R}}$, and it is straightforward to formulate a numerical algorithm of its solution. The system solved, the effective bending angle $ϵ$ equals the angle between ${\tilde{\mathbf{u}}}_{\mathrm{T}}$ and ${\tilde{\mathbf{u}}}_{\mathrm{R}}$ and the effective impact parameter $p=\left|{\mathbf{x}}_{\mathrm{T}}\times {\tilde{\mathbf{u}}}_{\mathrm{T}}\right|=\left|{\mathbf{x}}_{\mathrm{R}}\times {\tilde{\mathbf{u}}}_{\mathrm{R}}\right|$. It is important to realize that this is not an assumption or approximation. This is just the definition of the observables. Another way of expressing this is to say that the sphericity assumption errors in the observational data processing and the observation operator will mutually cancel out.

The physical model of RO observations is based on the GO ray Equation (6) complemented with the algorithm of the interpolation of gridded field of refractivity. However, this equation, by itself, allows the evaluation of a ray with the prescribed initial conditions. In a spherically layered atmosphere, this is sufficient to evaluate the impact parameter. However, in the presence of horizontal gradients, the impact parameter can only be evaluated from ${\tilde{\mathbf{u}}}_{\mathrm{T},\mathrm{R}}$, which are function of the relative Doppler frequency shift d, which is a function of ${\tilde{\mathbf{u}}}_{\mathrm{T},\mathrm{R}}$. One of the possibilities would be to solve the boundary problem, i.e., iterative solve for a ray connecting two satellites [60]. This may be conjugated with difficulties in the presence of multipath propagation, where there are multiple solutions.

We use a different approach. The observed bending angle profile ${ϵ}_{\mathrm{obs}}\left(p\right)$ satisfies the following geometrical equation:

$${\theta }_{RT}\left(t\right)={ϵ}_{\mathrm{obs}}\left(p\right)+\arccos \frac{p}{r_T\left(t\right)}+\arccos \frac{p}{r_R\left(t\right)},$$

(30)

where ${\theta }_{RT}\left(t\right)$ is the angle between vectors ${\mathbf{x}}_{\mathrm{T}}\left(t\right)$ and ${\mathbf{x}}_{\mathrm{R}}\left(t\right)$, and $r_{T,R}\left(t\right)=\left|{\mathbf{x}}_{\mathrm{T},\mathrm{R}}\left(t\right)\right|$. This equation can be numerically solved for the dependence of $t\left(p\right)$, which specifies the moment of time, when the ray with the impact parameter p was observed. The satellite positions for this moment of time will then define the occultation plane. In this plane we locate a ray starting at the transmitter and having the impact parameter p. The initial approximation is the starting ray direction ${\mathbf{u}}_{\mathrm{T}}$ satisfying the relation $p=\left|{\mathbf{x}}_{\mathrm{T}}\times {\mathbf{u}}_{\mathrm{T}}\right|$. In more detail we will discuss the iterative algorithm below, when discussing the adjoint version of the observation operator. Each ray is integrated to the point nearest to the receiver, because, generally speaking, the ray with effective impact parameter p and starting at the transmitter position ${\mathbf{x}}_{\mathrm{T}}$ will, most likely, not pass through the receiver position ${\mathbf{x}}_{\mathrm{R}}$, unless the model atmospheric state coincides with the actual atmospheric state.

3.1.2. Model of the 3-D Refractivity Field and Its Derivatives

To design a computational model of radio occultation experiments, a continuous model of 3-D refractivity field and its derivatives is required. In particular, the solution of the diffractive problem needs field $n\left(\mathbf{x}\right)$, the numerical integration of the geometric optical ray equation needs both $n\left(\mathbf{x}\right)$ and its gradient $\nabla n\left(\mathbf{x}\right)$, and the linear tangent model based on the perturbation theory requires Hessian matrix $\nabla \otimes \nabla n\left(\mathbf{x}\right)$. This can be done by interpolating the gridded field of the refractivity computed from the gridded fields of the specific model variables describing the atmospheric state in a specific model. For example, these variables may include temperature and humidity given at full and half sigma levels, and surface pressure, as in ECHAM models [14,61].

In a general case, there is a set of model variables $\mathbf{M}$ in the form of gridded fields of different model variables. For example, for a model with sigma levels, this vector will be represented as $\mathbf{M}=\left\{\left\{P_{jk}^s\right\},\left\{T_{ijk}\right\},\left\{q_{ijk}\right\}\right\}$, where $P_{jk}^s$ is the surface pressure at the surface grid with indexes $j,k$, $T_{ijk}$ is the gridded temperature, and $q_{ijk}$ is the gridded specific humidity with index i enumerating the vertical levels [61]. The surface grid may be a latitude-longitude grid ${\phi }_j,{\lambda }_k$ or more complicated, for example an icosahedral grid ${\phi }_{jk},{\lambda }_{jk}$ [62]. For the sake of simplicity, we will consider latitude-longitude grids.

For any model and any grid type, it is possible to evaluate the gridded field of refractivity $n_{ijk}=n(z_{ijk},{\phi }_j,{\lambda }_k)$ related to geometrical grids of altitudes $z_{ijk}$, latitudes ${\phi }_j$, and longitudes ${\lambda }_k$. This gridded field can then be interpolated to any spatial point. The interpolation procedure should be complemented with the extrapolation above the highest model grid, where we add mode altitude levels up to a height of about 120 km, and define the refractivity as a background model $n_{\mathrm{BG}}\left(z,\phi ,\lambda \right)$ multiplied with a fitting coefficient $\alpha \left(\phi ,\lambda \right)$, chosen to minimize the difference between the model temperature and pressure and those obtained by the hydrostatic integration of the merged refractivity profile.

Because refractive index $N=n-1$ decays with the height nearly exponentially, and for the ray integration smooth gradient of refractive index is necessary, we use the spline interpolation of $\ln N_{ijk}=\ln (n_{ijk}-1)$ as function of $z_{ijk}$ for each vertical profile, i.e., for each fixed pair of indexes $i,j$. For given point $\mathbf{x}$ in the Cartesian coordinates, its geodetic coordinates $(z,\phi ,\lambda )$ are calculated. Then horizontal grid mesh $({\phi }_J..{\phi }_{J+1},{\lambda }_K..{\lambda }_{K+1})$ containing point $(\phi ,\lambda )$ is located. Vertically interpolated values $\ln N_{jk}\left(z\right)$, $ln{}^{\prime }$$N{}_{jk}\left(z\right)$, $ln{}^{″}$$N{}_{jk}\left(z\right)$ for the four pairs of indexes $j=J..J+1$ and $k=K..K+1$ are calculated. The linear interpolation of these values with respect to $\phi ,\lambda$-coordinates is then performed to produce $\ln N\left(\mathbf{x}\right)$ and its derivatives $l{n}_z{}^{\prime }$$N\left(\mathbf{x}\right)$ and $l{n}_z{}^{″}$$N\left(\mathbf{x}\right)$ in the vertical direction.

Introducing the local vertical vector $\mathbf{v}\left(\mathbf{x}\right)=(\cos \phi \cos \lambda ,\cos \phi \sin \lambda ,\sin \phi )$ in the Cartesian coordinates, we can approximately calculate the gradient and the Hessian matrix:

$$\begin{array}{cc}\hfill \nabla n\left(\mathbf{x}\right)& =\mathbf{v}\left(\mathbf{x}\right)N\left(\mathbf{x}\right){\ln }_z^{\prime }N\left(\mathbf{x}\right),\hfill \\ \hfill \nabla \otimes \nabla n\left(\mathbf{x}\right)& =\mathbf{v}\left(\mathbf{x}\right)\otimes \mathbf{v}\left(\mathbf{x}\right)N\left(\mathbf{x}\right)\left[{\ln }_z^{″}N\left(\mathbf{x}\right)+{\left({\ln }_z^{\prime }N\left(\mathbf{x}\right)\right)}^2\right].\hfill \end{array}$$

(31)

The idea is to neglect the horizontal component of the gradient, which was found to be a good approximation. Taking into account the horizontal component of the gradient would require more complicated a horizontal interpolation, while the linear horizontal interpolation produces piecewise constant derivatives of $\ln N$ with respect to $\phi ,\lambda$, which are not continuous at the mesh borders.

3.1.3. Variations of Refractivity

The first part of the linear tangent, which is also used in the linear adjoint model, describes the variations of the refractivity due to variations of the model parameters. Since the ray trajectory equations only include the combination $n\nabla n$, which in our model is assumed to be equal to $\mathbf{v}\langle \mathbf{v},n\nabla n\rangle$, we only calculated the dependence of $\langle \mathbf{v},n\nabla n\rangle$ on the model variables, i.e., derivatives:

$$\frac{\partial \langle \mathbf{v}\left(\mathbf{x}\right),n\nabla n\left(\mathbf{x}\right)\rangle }{\partial \mathbf{M}}$$

(32)

The corresponding variations are then calculated as follows:

$$\delta \langle \mathbf{v}\left(\mathbf{x}\right),n\nabla n\left(\mathbf{x}\right)\rangle =\frac{\partial \langle \mathbf{v}\left(\mathbf{x}\right),n\nabla n\left(\mathbf{x}\right)\rangle }{\partial \mathbf{M}}\delta \mathbf{M}$$

(33)

The derivatives are evaluated by the differentiation of the interpolation scheme.

3.1.4. Variations of Ray Geometry

The geometric optical model is based on the numerical integration of ray trajectory Equation (6). Introducing augmented vector ${\displaystyle \mathbf{z}=\left(\genfrac{}{}{0pt}{}{\mathbf{x}}{\mathbf{u}}\right)}$ and denoting the right part of system (6) ${\displaystyle \mathbf{F}=\left(\genfrac{}{}{0pt}{}{\mathbf{u}}{n\nabla n\left(\mathbf{x}\right)}\right)}$, we rewrite the ray trajectory equation as follows:

$$\dot{\mathbf{z}}=\mathbf{F}\left(\mathbf{z}\right)$$

(34)

This equation is integrated numerically using finite steps. Remembering the definition of $\mathbf{F}\left(\mathbf{z}\right)$, we express its operator derivatives as follows:

$${\widehat{\mathbf{B}}}_m^{\mu }\equiv {\left.\frac{\partial \mathbf{F}\left(\mathbf{z}\right)}{\partial \mathbf{z}}\right|}_{\mathbf{z}={\mathbf{z}}_{m-1}^{\mu }}=\left(\begin{array}{cc}\widehat{\mathbf{0}}& \widehat{\mathbf{I}}\\ \nabla \otimes \nabla n^2\left({\mathbf{x}}_{m-1}^{\mu }\right)/2& \widehat{\mathbf{0}}\end{array}\right)\approx \left(\begin{array}{cc}\widehat{\mathbf{0}}& \widehat{\mathbf{I}}\\ \nabla \otimes \nabla n\left({\mathbf{x}}_{m-1}^{\mu }\right)& \widehat{\mathbf{0}}\end{array}\right)$$

(35)

where $\widehat{\mathbf{0}}$ and $\widehat{\mathbf{I}}$ are the zero and unit matrices of dimension $3\times 3$, index m enumerates integration steps, and index $\mu$ enumerates the substeps within each integration step, defined by the specific numerical integration scheme.

Introducing variations $\overline{\delta }{\mathbf{F}}_{m-1}^{\mu }$ of the form of the right part due to variations of the model refractivity, we can derive the following expression for variations of ${\mathbf{z}}_m$:

$$\delta {\mathbf{z}}_m={\widehat{\mathbf{B}}}_m\delta {\mathbf{z}}_{m-1}+\sum _{\mu }{\widehat{\mathbf{C}}}_m^{\mu }\overline{\delta }{\mathbf{F}}_{m-1}^{\mu }$$

(36)

where matrices ${\widehat{\mathbf{B}}}_m$ and ${\widehat{\mathbf{C}}}_m^{\mu }$ are obtained by differentiating the specific numerical integration scheme and are expressed as polynomial functions of matrices ${\widehat{\mathbf{B}}}_m^{\mu }$. An example of the complete evaluation of these matrices for the Runge–Kutta scheme of the fifth order can be found in [14].

Introducing notation ${\alpha }_{m-1}^{\mu }=\langle \mathbf{v}\left({\mathbf{x}}_{m-1}^{\mu }\right),n\nabla n\left({\mathbf{x}}_{m-1}^{\mu }\right)\rangle$ for the parameters influencing the ray geometry, and uniting parameters ${\alpha }_{m-1}^{\mu }$ into vector ${\mathbf{a}}_{m-1}$, we arrive at the expression for the variations of the form of the right part:

$$\overline{\delta }{\mathbf{F}}_{m-1}^{\mu }={\widehat{\mathbf{A}}}_m^{\mu }\delta {\mathbf{a}}_{m-1}={\widehat{\mathbf{A}}}_m^{\mu }\frac{{\partial \mathbf{a}}_{m-1}}{\partial \mathbf{M}}\delta \mathbf{M}$$

(37)

where the matrices ${\widehat{\mathbf{A}}}_m^{\mu }$ cut out ${\alpha }_{m-1}^{\mu }$ from ${\mathbf{a}}_{m-1}$. Collecting all these transformations, we arrive at the equation for the variations of augmented vector $\mathbf{z}$:

$$\delta {\mathbf{z}}_m={\widehat{\mathbf{B}}}_m\delta {\mathbf{z}}_{m-1}+\sum _{\mu }{\widehat{\mathbf{C}}}_m^{\mu }{\widehat{\mathbf{A}}}_m^{\mu }\frac{{\partial \mathbf{a}}_{m-1}}{\partial \mathbf{M}}\delta \mathbf{M}\equiv {\widehat{\mathbf{B}}}_m\delta {\mathbf{z}}_{m-1}+{\widehat{\mathbf{C}}}_m\delta \mathbf{M}$$

(38)

3.1.5. Variations of Refraction Angle

The refraction angle and impact parameter are functions of the initial and final conditions of a ray trajectory, as defined by Equation (29):

$$\begin{array}{c}\hfill ϵ=ϵ({\mathbf{z}}_0,{\mathbf{z}}_N)\\ \hfill p=p({\mathbf{z}}_0,{\mathbf{z}}_N)\end{array}$$

(39)

Final condition ${\mathbf{z}}_N$ is a function of the initial condition and model variables:

$${\mathbf{z}}_N={\mathbf{z}}_N({\mathbf{z}}_0,\mathbf{M})$$

(40)

The full variations of the refraction angle and impact parameter can be written in the following form:

$$\begin{array}{c}\hfill \delta ϵ=\frac{\partial ϵ}{\partial {\mathbf{z}}_0}\delta {\mathbf{z}}_0+\frac{\partial ϵ}{\partial {\mathbf{z}}_N}\frac{\partial {\mathbf{z}}_N}{\partial {\mathbf{z}}_0}\delta {\mathbf{z}}_0+\frac{\partial ϵ}{\partial {\mathbf{z}}_N}\frac{\partial {\mathbf{z}}_N}{\partial \mathbf{M}}\delta \mathbf{M}\\ \hfill \delta p=\frac{\partial p}{\partial {\mathbf{z}}_0}\delta {\mathbf{z}}_0+\frac{\partial p}{\partial {\mathbf{z}}_N}\frac{\partial {\mathbf{z}}_N}{\partial {\mathbf{z}}_0}\delta {\mathbf{z}}_0+\frac{\partial p}{\partial {\mathbf{z}}_N}\frac{\partial {\mathbf{z}}_N}{\partial \mathbf{M}}\delta \mathbf{M}\end{array}$$

(41)

We need the variation of the refraction angle with a given impact parameter, so we choose the variation of the initial condition so that $\delta p$ should be equal to 0. To arrive at a completely determined system of conditions for $\delta {\mathbf{z}}_0$, we assume that we only vary ray direction ${\mathbf{u}}_0$, and that its variation is coplanar with vectors ${\mathbf{x}}_0$ and ${\mathbf{x}}_N$. Complementing this with the requirement for varied ${\mathbf{u}}_0$ to remain a unit vector, we can uniquely define $\delta {\mathbf{z}}_0$ from the following system:

$$\begin{array}{c}{\displaystyle \left({\displaystyle \frac{\partial p}{\partial {\mathbf{z}}_0}+\frac{\partial p}{\partial {\mathbf{z}}_N}\frac{\partial {\mathbf{z}}_N}{\partial {\mathbf{z}}_0}}\right)\delta {\mathbf{z}}_0\equiv \frac{dp}{d{\mathbf{z}}_0}\delta {\mathbf{z}}_0=-\frac{\partial p}{\partial {\mathbf{z}}_N}\frac{\partial {\mathbf{z}}_N}{\partial \mathbf{M}}\delta \mathbf{M}}\hfill \\ \delta {\mathbf{x}}_0=0\hfill \\ (\delta {\mathbf{u}}_0,[{\mathbf{x}}_0,{\mathbf{x}}_N])=0\hfill \\ ({\mathbf{u}}_0,\delta {\mathbf{u}}_0)=0\hfill \end{array}$$

(42)

The full derivative ${\displaystyle \frac{dp}{d{\mathbf{z}}_0}}$ is used in the iterative solution for the initial condition of ray with prescribed impact parameter p. We shall use the following symbolic notation for the solution of this system:

$$\delta {\mathbf{z}}_0=-\frac{d{\mathbf{z}}_0}{dp}\frac{\partial p}{\partial {\mathbf{z}}_N}\frac{\partial {\mathbf{z}}_N}{\partial \mathbf{M}}\delta \mathbf{M}$$

(43)

For the variation of the refraction angle, we have the following expression:

$$\delta ϵ=\left[\frac{\partial ϵ}{\partial {\mathbf{z}}_N}-\left(\frac{\partial ϵ}{\partial {\mathbf{z}}_0}+\frac{\partial ϵ}{\partial {\mathbf{z}}_N}\frac{\partial {\mathbf{z}}_N}{\partial {\mathbf{z}}_0}\right)\frac{d{\mathbf{z}}_0}{dp}\frac{\partial p}{\partial {\mathbf{z}}_N}\right]\frac{\partial {\mathbf{z}}_N}{\partial \mathbf{a}}\delta \mathbf{a}\equiv \frac{dϵ}{d{\mathbf{z}}_N}\frac{\partial {\mathbf{z}}_N}{\partial \mathbf{M}}\delta \mathbf{M}$$

(44)

The symbolic full derivative ${\displaystyle \frac{dϵ}{d{\mathbf{z}}_N}}$ introduced here describes the sensitivity of the refraction angle with respect to the ray geometry. Equation (44) is the expression of the adjoint observation operator for the bending angle as a function of the impact parameter.

3.1.6. Error Covariances

The assimilation requires estimates of error covariances. The ionospheric correction algorithm combined with the statistical optimization and with the residual error estimate was developed by [51]. A simple model of covariances is described by [63]. A dynamical algorithm of the covariance estimate was introduced by [64]. Radio holographic estimates of bending angles errors were described by [65] and further developed by [66]. More advanced analysis of uncertainty propagation through the retrieval chain and resulting covariances was performed by [53,54,67,68,69].

3.2. 1-D Assimilation of Bending Angles

3.2.1. Operators Based on the Abel Integral

The 1-D assimilation of bending angles [17,18,19,20,21,22,23] uses the representation of the bending angle as the Abel integral of a single profile of refractivity gradient over the refractive radius $x=nr$, as specified by Equation (23). This representation is convenient, because $ϵ\left(p\right)$ is a linear functional of $n\left(x\right)$.

There are different numerical schemes of evaluation of the Abel integral on finite grids. Ref. [20] use a piecewise exponential representation of $N\left(x\right)$ between the nodes of grid $x_i,i=1,...,K$:

$$\begin{array}{cc}\hfill & \ln n\approx n-1\equiv N,\hfill \\ \hfill & k_i=\frac{\ln \left(N_i/N_{i+1}\right)}{x_{i+1}-x_i},\hfill \\ \hfill & \frac{d\ln n}{dx}\approx k_i{N}_i\exp \left(-k_i\left(x-x_i\right)\right),x_i\le x\le x_{i+},\hfill \end{array}$$

(45)

as well as the following approximation of the integral kernel:

$$\sqrt{x^2-p^2}\approx \sqrt{2p}\sqrt{x-p}.$$

(46)

This results in the following discrete approximation of the bending angle:

$$ϵ\left(p\right)=\sqrt{2\pi p}\sum _{i=1}^{K-1}{k}_i{N}_i\exp \left(k_i\left(x_i-p\right)\right)\left[\mathrm{erf}\left(\sqrt{k_i\left(x_{i+1}-p\right)}\right)-\mathrm{erf}\left(\sqrt{k_i\left(x_i-p\right)}\right)\right].$$

(47)

Alternatively, the following coordinate change [21] is performed:

$$x=\sqrt{p^2+s^2},$$

(48)

which eliminates the pole of the integral kernel and allows the numerical integration in an equally spaced grid in s, using the trapezoidal rule.

The linear adjoint model uses the block for the evaluation of the derivatives of the refractivity gradient with respect to the model variables. It must be complemented with the following relationships:

$$\begin{array}{cc}\hfill & \frac{dn}{dx}=\frac{{n^{\prime }}_r}{n+r{n^{\prime }}_r},\hfill \\ \hfill & {\left.\frac{\partial }{\partial \mathbf{M}}\frac{dn}{dx}\right|}_x=\frac{\partial }{\partial \mathbf{M}}\frac{dn}{dx}-\frac{d^{2}n}{d{x}^2}\frac{\partial x}{\partial \mathbf{M}},\hfill \\ \hfill & {\left.\frac{\partial n}{\partial \mathbf{M}}\right|}_x=\frac{\partial n}{\partial \mathbf{M}}-\frac{dn}{dx}\frac{\partial x}{\partial \mathbf{M}},\hfill \end{array}$$

(49)

the two last expressions providing the derivatives for a fixed value of x. They are necessary, because $x=nr$ also depends on the refractivity.

These approaches based on the Abel integral over refractive radius x rely upon $n\left(x\right)$ being a single-valued function. The above discussion of the waveguides in Section 2.3 indicates that for the points located below a waveguide, bending angle, although being defined, requires the integration over two branches of a multi-valued dependence $n\left(x\right)$ or the integration over the geometrical radius r, expressed by Equation (22). In this case, however, simple analytical formulas for the bending angle increment at each step of the grid cannot be derived due to a non-linear dependence of the sub-integral expression from n.

3.2.2. 1-D Ray-Tracing Operators

A straightforward solution for 1-D observation operator that can be applied in the presence of waveguides is the use of ray-tracing. This approach mimics the 3-D ray-tracing approach, as described in Section 3.1, with one single modification. Instead of using 3-D gridded fields of refractivity, one single 1-D vertical profile without horizontal interpolation is used. This approach will be referred to as 3D_CAR_RT, which stays for 3-D ray-tracing with constant atmospheric refractivity along horizontal. The 3-D ray-tracing taking into account the horizontal gradients of refractivity will be referred to as 3D_RT. From the viewpoint of the computational expenses, 3D_CAR_RT approach is much cheaper than 3D_RT. On the other hand, 3D_CAR_RT is not equivalent to the Abel integral, because it can simulate the realistic geometry of the Earth, including the reference ellipsoid and geoid, rather than spherical symmetry centered at the local curvature center.

Instead of using the Abel integral, it is possible to use the special form of the ray trajectory equation for spherically layered medium, from which the Abel integral follows, as shown in Section 2. We will refer to this approach as 1D_RT. The simplest equation with the order reduced by using a symmetry is Equation (14) written for a single half of the ray trajectory:

$$\frac{dr}{d\theta }=\frac{r}{p}\sqrt{n^2{r}^2-p^2}.$$

(50)

This equation does not contain the gradient of refractivity, which makes it numerically more stable. Due to this, 1D_RT has a smaller computational demand as compared to 3D_CAR_RT, although both these approaches integrate a differential equation depending on just one vertical profile of refractivity.

This equation, however, needs some special treatment of the perigee point: given the perigee radius $r_p$ that satisfies the equation $n\left(r_p\right)r_p=p$, Equation (50) has a constant solution $r\left(\theta \right)=r_p$. Therefore, it must be complemented with the initial condition for the second derivative. Using Equation (15) and relation $\Delta \theta =n\Delta t/r$ at the perigee point, we start the integration at $r\left(\theta =0\right)=r_p$, and make the first step:

$$r\left(\Delta \theta \right)=r_p+\frac{{\Delta \theta }^2}{2}{\left.\frac{r}{n}\frac{dx}{dr}\right|}_{r=r_p}=r_p+\frac{{\Delta \theta }^2}{2}{r}_p\left(1+r_p\frac{n^{\prime }}{n}\right).$$

(51)

After that, Equation (50) is integrated, until r reaches the pre-specified upper boundary, e.g., $r_E+120$ km, where $r_E$ is the Earth’s local curvature radius. The bending angle is then evaluated as follows:

$$ϵ=2\theta +2\arcsin \frac{p}{r\left(\theta \right)}-\pi ,$$

(52)

assuming that $n\left(r\right)=1$. The tangent linear operator follows the same guidelines as in Section 3.1, but is much simpler. First, the ray is integrated for a prescribed value of impact parameter p. Instead of 6 dynamic variables, there is only one. Differentiating the numerical integration scheme of Equation (50), with the initial condition:

$$\frac{\partial r_p}{\partial \mathbf{M}}=-\frac{\partial n}{\partial \mathbf{M}}\frac{r_p}{n+r_p{n}_r^{\prime }},$$

(53)

similarly to Equation (44), we write

$$\delta ϵ=\frac{dϵ}{dr}\frac{\partial r\left(\theta \right)}{\partial \mathbf{M}}\delta \mathbf{M}=\frac{p}{r\left(\theta \right)\sqrt{r{\left(\theta \right)}^2-p^2}}\frac{\partial r\left(\theta \right)}{\partial \mathbf{M}}\delta \mathbf{M}.$$

(54)

Here the dependence on the model state vector $\mathbf{M}$ only includes one vertical profile.

3.3. Assimilation of Refractivity

3.3.1. 1-D Assimilation of Refractivity

1-D assimilation of refractivity practically does not differ from 1-D assimilation of bending angles, because refractivity profile is retrieved using the sphericity assumption and, therefore, is uniquely linked to the bending angle profile. It is only necessary to evaluate the covariances for the refractivity [70].

3.3.2. Refractivity-Mapping Operators

Refractivity-mapping operators [25,26,27,28,29,30] should be looked at as the 3-D assimilation of refractivity. The basic idea of refractivity mapping is to evaluate the 1-D retrieved refractivity profile as a functional of 3-D refractivity field. Using an exact form of this functional would not bring anything new as compared to the 3-D assimilation of bending angles, in a way similar to 1-D assimilation of refractivity. The further idea is to use some approximate solutions for the ray trajectories that make this operator less computationally expensive.

The idea originates from [71], who derived the 2-D resolution kernel. The starting point of [71] was the approach used by [72] and originally developed by [73]. Equation (20) allows expressing the bending angle as an integral of the refractivity gradient multiplied with the ray normal vector over the ray. The exact evaluation of this integral requires the integration of the ray trajectory equations. However, using the straight-line approximation [25], neglecting the horizontal component of the gradient, and approximating logarithm of refractivity $\ln n$ as refractive index $n-1=N$ it can be written as follows:

$$ϵ\left(p\right)=p\underset{p}{\overset{\infty }{\int }}\frac{\partial N\left({\displaystyle r,\pm \arccos \frac{p}{r}}\right)}{\partial r}\frac{dr}{\sqrt{r^2-p^2}},$$

(55)

where $N=N\left(r,\theta \right)$ is a function of the radial coordinate r and the angular coordinate $\theta$ in the occultation plane, $\theta =0$ corresponding to the ray perigees. This expression is understood as the sum over two parts of the ray trajectory, which accounts for the missing factor of 2 in front of the integral. Hereinafter, expressions containing ± and/or ∓ signs are also understood in this sense. The corresponding Abel inversion has the following form:

$$\tilde{N}\left(r\right)=-\frac{1}{\pi }\underset{r}{\overset{\infty }{\int }}\frac{ϵ\left(p\right)dp}{\sqrt{p^2-r^2}}.$$

(56)

Substituting Equation (55) into the Equation (56), we arrive at the following operator composition:

$$\tilde{N}\left(r\right)=-\frac{1}{\pi }\underset{r}{\overset{\infty }{\int }}\underset{r}{\overset{r^{\prime }}{\int }}\frac{\partial N\left({\displaystyle r^{\prime },\pm \arccos \frac{p}{r^{\prime }}}\right)}{\partial r^{\prime }}\frac{pdp}{\sqrt{p^2-r^2}\sqrt{{r^{\prime }}^2-p^2}}d{r}^{\prime }.$$

(57)

where $\tilde{N}\left(r\right)$ is the retrieved refractivity. If $N=N\left(r\right)$, we can evaluate the integral over p using the following identity:

$$\underset{r}{\overset{r^{\prime }}{\int }}\frac{pdp}{\sqrt{p^2-r^2}\sqrt{{r^{\prime }}^2-p^2}}=\frac{\pi }{2},$$

(58)

and arrive at the expected expression for a spherically layered medium:

$$\tilde{N}\left(r\right)=-\underset{r}{\overset{\infty }{\int }}\frac{dN\left(r^{\prime }\right)}{d{r}^{\prime }}d{r}^{\prime }=N\left(r\right).$$

(59)

Changing in Equation (57) the coordinate p to $\theta =\arccos \left(p/r^{\prime }\right)$, we can rewrite it in a simple form:

$$\tilde{N}\left(r\right)=-\frac{1}{\pi }\underset{r}{\overset{\infty }{\int }}\underset{0}{\overset{\arccos \frac{r}{r^{\prime }}}{\int }}\frac{\partial N\left(r^{\prime },\pm \theta \right)}{\partial r^{\prime }}\frac{\cos \theta d\theta }{\sqrt{{\cos }^2\theta -{\left({\displaystyle \frac{r}{r^{\prime }}}\right)}^2}}d{r}^{\prime }.$$

(60)

To make further simplifications, we represent $N\left(r^{\prime },\theta \right)$ as follows:

$$N\left(r,\theta \right)=N\left(r,0\right)+\underset{0}{\overset{\theta }{\int }}\frac{\partial N\left(r,{\theta }^{\prime }\right)}{\partial {\theta }^{\prime }}d{\theta }^{\prime }.$$

(61)

Substituting this into Equation (60) and integrating by parts (or switching the integration order), we arrive at the following expression:

$$\tilde{N}\left(r\right)=N\left(r\right)\mp \frac{1}{\pi }\underset{r}{\overset{\infty }{\int }}\underset{0}{\overset{\arccos \frac{r}{r^{\prime }}}{\int }}\frac{{\partial }^{2}N\left(r^{\prime },\pm {\theta }^{\prime }\right)}{\partial r^{\prime }\partial {\theta }^{\prime }}\underset{{\theta }^{\prime }}{\overset{\arccos \frac{r}{r^{\prime }}}{\int }}\frac{\cos \theta d\theta }{\sqrt{{\cos }^2\theta -{\left({\displaystyle \frac{r}{r^{\prime }}}\right)}^2}}d{\theta }^{\prime }d{r}^{\prime }.$$

(62)

Analytically evaluating the integral over $\theta$ and changing notation from ${\theta }^{\prime }$ to $\theta$, we finally arrive at the following expression:

$$\begin{array}{cc}\hfill \tilde{N}\left(r\right)& =N\left(r,0\right)\mp \frac{1}{\pi }\underset{r}{\overset{\infty }{\int }}\underset{0}{\overset{\arccos \frac{r}{r^{\prime }}}{\int }}\frac{{\partial }^{2}N\left(r^{\prime },\pm \theta \right)}{\partial r^{\prime }\partial \theta }\arccos \frac{r^{\prime }\sin \theta }{\sqrt{{r^{\prime }}^2-r^2}}d\theta d{r}^{\prime }\hfill \\ \hfill & =N\left(r,0\right)\mp \frac{1}{\pi }\underset{0}{\overset{\pi /2}{\int }}\underset{r/\cos \theta }{\overset{\infty }{\int }}\frac{{\partial }^{2}N\left(r^{\prime },\pm \theta \right)}{\partial r^{\prime }\partial \theta }\arccos \frac{r^{\prime }\sin \theta }{\sqrt{{r^{\prime }}^2-r^2}}d{r}^{\prime }d\theta .\hfill \end{array}$$

(63)

This expression represents the mapping operator as the sum of the local refractivity profile and correction term depending on the horizontal gradient. This representation of resolving kernel is simpler than that derived by [71]. Although this expression was derived using the straight-line approximation, the approximation errors will only apply to non-spherical component of the refractivity field, because the spherically symmetrical component will transform identically. On the other hand, Ref [25] suggested a simple way of taking realistic ray trajectories into account. Ray trajectories can be evaluated for the first-guess refractivity field, and then new local coordinate ${\theta }^{\prime }$ in the occultation plane can be introduced, which takes into account the ray deviation from the straight-line ray equation.

4. Numerical Results

In this section, we compare the performances of the three different RO operators using Spire measurements collected between 8–28 March 2018. Spire data is generated by a constellation of low Earth orbit CubeSats in contrast with the more traditional larger satellites such as METOP and COSMIC. Bias and standard deviation of relative bending angles errors are computed using 5377 Spire RO profiles. Spire data employs an equidistant vertical grid using increments of 0.2 km. Because at altitudes higher than 50 km, the profiles are seriously affected by ionospheric errors, we decided to limit the comparison up to this altitude threshold. Finally, the statistics are computed using 4 cycles per day and 6 h GFS forecasts as background.

Figure 2 and Figure 3 show the bias and standard deviation for Spire BA data, evaluated for three different observation operators, denoted as follows: 1D_RT is 1-D ray-tracing based on (50), 3D_CAR_RT is the 3D ray-tracing with constant atmospheric refractivity along horizontal, and 3D_RT is the full 3-D ray-tracing.

5. Discussion

The full 3-D ray tracing observation operator, in most cases, reduces the systematic differences between GFS forecasts and RO observations in the lower troposphere, as compared to the other two simplified observation operators. The larger systematic difference for the simplest 1-D ray tracing operator at large heights come from numerical integration errors and can be reduced by using a smaller integration step. In the lower troposphere, especially in the tropics, the lowest standard deviation is achieved for the simplest 1-D ray tracing, while for the full 3-D ray tracing it is the highest. To better understand this, a further impact study is required.

From the viewpoint of the computational expenses, 3D_CAR_RT approach is much cheaper than 3D_RT. This follows from the fact tha the highest computational expenses come from the evaluation of its derivatives with respect to the gridded fields of the model variables. Therefore, in 3D_RT, the evaluation of $\partial n/\partial \mathbf{M}$ involves much more components of the atmospheric state vector $\mathbf{M}$, as compared to 3D_CAR_RT. On the other hand, 3D_RT approach can be limited to the lower troposphere and combined with another 1-D operator above a height of 7–10 km, where the influence of the horizontal gradients upon the standard retrieval based on the Abel inversion becomes weaker.

3D_CAR_RT and 1D_RT approaches are very much equivalent in terms of computational expenses. In both cases, 1-D integration is used, and the evaluation of $\partial n/\partial \mathbf{M}$ only involve a single vertical profile of model state. Still, 3D_CAR_RT is not equivalent to the Abel integral, because it is capable of simulating the realistic geometry of the Earth, including the reference ellipsoid and geoid, rather than spherical symmetry centered at the local curvature center.

The refractivity-mapping operator (63) can be looked at as a 3-D solution that optimizes the computational costs. Although it requires the knowledge of derivatives $\partial n/\partial \mathbf{M}$ taken with respect to the full 3-D atmospheric state vector, they only need to be evaluated at a regular spatial grid, unlike the 3-D raytracing, where the derivatives are evaluated along different rays. Therefore, using an approximation of fixed occultation plane, it is possible to evaluate the necessary set of derivatives at a regular grid just once.

6. Conclusions

In this paper, we discussed different strategies of RO data assimilation into models of global atmospheric circulation. As was early recognized, the most convenient variable to be assimilated is the bending angle. Because the application of data processing techniques based on Fourier Integral Operators, implementing canonical transforms in the wave optics, allows the retrieval of GO bending angles, the problem can be reduced to the geometrical optics. It should also be possible to assimilate the retrieved refractivity, but it is equivalent to the assimilation of bending angles, because the refractivity is uniquely linked to the bending angle by the Abel transform. The basic equations for the bending angle can be directly derived from the wave equation, which can be solved under the GO approximation. The GO ray equations in different forms are the basis for the formulation of observation operators. When processing real RO observation, we cannot derive the bending angle and ray impact parameter, because their derivation requires an additional assumption of the spherical symmetry. Because the real atmosphere is not spherically symmetric, we introduce the concepts of the effective impact parameter and bending angle, which are defined in the way they are evaluated in the spherically symmetric case. This definition is then included into the observation operator, and, therefore, any errors due the non-sphericity of the atmosphere cancel out. We discuss different approaches to the formulation of the BA observation operator. The most efficient numerical implementation is based on the 1-D scheme, where a single atmospheric profile is used and the atmosphere is assumed to be spherically symmetric with respect to the local curvature radius of the reference ellipsoid. The observation operator can be based on the Abel integral, but in this case it is susceptible to errors due to wave guides. We introduce an alternative formulation based on a reduced-order ray equation, from which the Abel integral can be derived. Still, this form of observation operator can correctly take waveguides into account, because it does not use the assumption that the refractive radius is uniquely linked to the geometric radius. We discuss the most accurate observation operator based on the full 3-D ray-tracing equations. This operator is computationally expensive. However, it is possible to optimize the computation power demand by only using this operator in the lower troposphere, for heights below 7 km. The 3-D ray-tracing can be in a straightforward way reformulated as 1-D ray-tracing model. To this end, instead of using a 3-D interpolation scheme of the gridded refractivity field, a 1-D interpolation that only involves a single vertical profile, is employed. The form of the observation operator based on the 3-D ray-tracing with constant atmospheric refractivity along the horizontal is more accurate than the 1-D operator based on the reduced form of ray trajectories in the spherically symmetric medium, because in the 3-D operator the realistic shape of the Earth is taken into account. We derived a simple analytic solution for the refractivity-mapping operator. This operator combines the numerical efficiency and the account of the horizontal gradients. We tested the three different variants of ray-tracing observation operators using the Spire CubeSat observations. The observation operator was evaluated for GFS forecast fields. It was shown that the full 3-D ray-tracing observation operator reduces the systematic differences between GFS forecasts and RO observations in the lower troposphere, as compared to the other two simplified observation operators. The larger systematic difference for the simplest 1-D ray-tracing operator occurs at large heights of Spire data, in the lower troposphere, especially in the tropics, the lowest standard deviation is achieved for the simplest 1-D ray-tracing, while for the full 3-D ray-tracing it is the highest.