3-D tomographic limb sounder retrieval techniques: irregular grids and Laplacian regularisation

2.1 The retrieval

The main focus of this paper is the algorithm for retrieving atmospheric quantities, such as temperature or trace gas mixing ratios, from limb sounder measurements in 3-D. In this section, we outline the general inverse modelling approach to this problem and identify some aspects we aim to improve upon.

Let an atmospheric state vector x∈ℝⁿ represent the values of some atmospheric quantities on a finite grid and let y∈ℝ^m be a set of m remote measurements taken within the region in question. The physics of the measurement process itself has to be understood for the instrument to be practical, so one can usually build a theoretical model of the measurement (the forward model) $F:{\mathbb{R}}^n\to {\mathbb{R}}^m$. The inverse problem of determining x given y is then typically both underdetermined and ill posed: many atmospheric states would result in the same measurement, but due to instrument and model errors no states would yield the exact measurement result y. One can solve this problem by finding an atmospheric state x that minimises the following quadratic cost function:

The first term of Eq. (1) quantifies the difference between the actual measurement and a simulated one (given the atmospheric state is x). The matrix ${\mathbf{S}}_{\mathit{ϵ}}^{-\mathrm{1}}\in {\mathbb{R}}^m\times {\mathbb{R}}^m$ represents the expected instrument errors. The second term introduces regularisation; i.e. it is larger for x that are unlikely given our general, measurement-unrelated knowledge about the atmosphere. In practice this is achieved using an a priori state x_a, usually derived from climatologies. It is meant for introducing the best of our prior knowledge into the retrieval process without imposing any features of the kind to be measured. The precision matrix ${\mathbf{S}}_{\mathrm{a}}^{-\mathrm{1}}$ (the inverse of a covariance matrix for atmospheric states S_a) contains prior knowledge about the probability distribution of atmospheric states. For a more detailed discussion of this approach refer to Rodgers (2000).

Finding the precision matrix ${\mathbf{S}}_{\mathrm{a}}^{-\mathrm{1}}$ that would correctly represent the physical and statistical properties of the atmosphere is challenging. The Tikhonov regularisation approach (Tikhonov and Arsenin, 1977)

is widely used for this purpose. Here the matrix L₀ is diagonal and contains reciprocals of standard deviations of atmospheric quantities, and L_x, L_y, L_z represent spatial derivatives in the respective directions on a regular rectangular grid, estimated by forward differences. The positive constants α₀, α_h, α_v are chosen ad-hoc. This is a convenient way to construct ${\mathbf{S}}_{\mathrm{a}}^{-\mathrm{1}}$, and it serves the purpose of ruling out very pathological, oscillatory retrievals. The constants α₀, α_h, α_v are, however, unphysical, grid dependent, and hence have to be estimated (typically by trial and error) for each use case. They are not related to the properties of the atmosphere in any simple fashion. Also, the usage of first-order derivatives only is justified mostly by simplicity and the need to keep the computations cheap.

We chose a slightly different, more physically and statistically motivated approach to regularisation, which is introduced in the following section.

2.2 Covariance

This section describes the theoretical background for the regularisation algorithm used for our retrieval, i.e. the motivation for ${\mathbf{S}}_{\mathrm{a}}^{-\mathrm{1}}$ in Eq. (1).

If the atmospheric state x has the mean x_apr and covariance matrix S_a then its entropy is maximised if x has multivariate normal distribution (proved e.g. by Rodgers, 2000). Hence, by treating x as a random vector with this distribution we impose the smallest restriction possible while still introducing a priori information from S_a into our retrieval. Also, the probability density of x is then

This immediately justifies the second term of Eq. (1), provided that ${\mathbf{S}}_{\mathrm{a}}^{-\mathrm{1}}$ represents, to some extent, the actual statistics of the atmosphere. This is often referred to as the Bayesian approach to regularisation (Rodgers, 2000). The precision matrix from Eq. (2) is, however, a mathematical device not meant to represent the physical world, as discussed in the previous section. Furthermore, we would like to derive a coordinate independent expression for the regularisation term of cost function (also unlike Eq. 2), as this would be useful for work with irregular grids and allow us to use the same regularisation on completely different grid geometries. To achieve these goals, we use continuous covariance operators and their associated norms and only discretise them as a last step, hence preserving the grid independence and their statistical interpretation. The resulting discrete regularisation is a variant of Tikhonov in its final numerical form, but can be interpreted as a realisation of general continuous covariance relations. The following paragraph introduces some of the required formalism (Lim and Teo, 2009; Tarantola, 2013).

Let us denote the departure of the atmospheric quantity f(r) from the a priori by $\mathit{\varphi }\left(\mathit{r}\right)=f\left(\mathit{r}\right)-f_{\mathrm{apr}}\left(\mathit{r}\right)$ for some r∈ℝ³. In some finite volume V⊂ℝ³, we define the covariance operator C by

where $C_{\mathrm{k}}:{\mathbb{R}}^{\mathrm{3}}\times {\mathbb{R}}^{\mathrm{3}}\to \mathbb{R}$ is the covariance kernel (also known as covariance function). Then we can treat the scalar fields ϕ(r) as elements of Hilbert space with the product

which induces the norm $\Vert \mathit{\varphi }{\Vert }^{\mathrm{2}}=\langle \mathit{\varphi },\mathit{\varphi }\rangle$. Once the explicit expression for the norm is found it can be discretised, representing the scalar field ϕ(r) by the atmospheric state vector x−x_apr. Then we can write

It now remains to find an appropriate covariance kernel C_k. Let us first consider an atmosphere that is isotropic and has the same physical and statistical properties everywhere. Then one would expect $C_{\mathrm{k}}\left(\mathit{r},{\mathit{r}}^{\mathbf{\prime }}\right)=g\left(\Vert \mathit{r}-{\mathit{r}}^{\mathbf{\prime }}\Vert \right)$ with some monotonously decreasing function $g:[\mathrm{0},\mathrm{\infty })\to [\mathrm{0},\mathrm{\infty })$. The most common kernel used in literature for fluid dynamics, meteorology and similar applications is the parametric Matérn covariance kernel (introduced, e.g. by Lim and Teo, 2009):

where $d=∥\mathit{r}-{\mathit{r}}^{\prime }∥$, σ and L are the standard deviation and typical length scale of structures (correlation length), respectively, of the atmospheric quantity in question. K_ν is the modified Bessel function of the second kind, Γ(ν) is the Gamma function. Exponential and Gaussian covariance kernels are special cases of the kernel (Eq. 7), with ν=0.5 and ν=1 respectively. We choose the exponential covariance

as it closely resembles the Matérn with ν values that are typically used for fluid problems and allows for analytic derivation of the subsequently required quantities. Also, having a significant number of parameters that could only be estimated theoretically, we could not make good use of the flexibility provided by the free parameter of the Matérn covariance in any case.

It can be shown (Tarantola, 2013) that the norm associated with the covariance (Eq. 8) is

neglecting some boundary terms. In a more realistic picture of the atmosphere, the correlation length L depends on altitude and strongly depends on direction: the correlation between vertically separated air parcels is much weaker than between air parcels separated by the same distance horizontally. We propose to deal with anisotropic or variable L by performing a coordinate transformation such that L would be isotropic and constant in resulting coordinates. In particular, let $U,V\subset {\mathbb{R}}^{\mathrm{3}}$ and consider a bijective map $\mathit{\xi }:V\to U$ such that both ξ and ξ⁻¹ are twice differentiable on their respective domains and have non-zero first and second derivatives everywhere. Then, using integration by substitution and basic vector calculus identities, Eq. (9) can be written as

where

Here $(\mathbf{D}\mathit{\xi }{)}_{ij}=\partial {\mathit{\xi }}_i/\partial u_j$ is the Jacobian, and we define the matrix $({\mathbf{D}}^{\mathrm{2}}f{)}_{ij}=({\partial }^{\mathrm{2}}f)/(\partial u_i\partial u_j)$. ξ can be chosen so that all its spatial derivatives could be computed analytically, and ∇(ϕ(ξ)), D²(ϕ(ξ)) are the derivatives of ϕ in the transformed space, so numerical evaluation of Eq. (10) is similar in complexity to that of Eq. (9), with the only notable difference being the need to compute the mixed derivatives ${\partial }^{\mathrm{2}}\mathit{\varphi }/(\partial u_i\partial u_j),i\ne j$ if Dξ is not diagonal. A large and relevant class of suitable maps ξ can be defined, for example, by assuming that the correlation lengths in vertical and horizontal direction are smooth functions of altitude L_v(z), L_h(z) (we use Cartesian coordinates $(x,y,z)$ where the z axis is vertical). Then the map

defines a corresponding transformation. If, in addition, L_h is constant, mixed derivatives are not required and we can simply evaluate ∇(ϕ(ξ)), Δ(ϕ(ξ)) in transformed space. Transformations of this type could be used, for example, to perform regularisation in tropopause-based coordinates and use different correlation lengths for troposphere and stratosphere. For the study described in this paper, we only introduce anisotropy (i.e. L_v and L_h are different, but do not depend on altitude) and Eq. (10) becomes

This general expression can now be discretised for use in the retrieval. We can represent the scalar field ϕ(r) by the atmospheric state vector x−x_apr and then calculate a discrete approximation to the integral in Eq. (13). This can then be interpreted as ${\left(\mathit{x}-{\mathit{x}}_{\mathrm{apr}}\right)}^T{\mathbf{S}}_{\mathrm{a}}^{-\mathrm{1}}\left(\mathit{x}-{\mathit{x}}_{\mathrm{apr}}\right)$ as seen from Eq. (6). For example, if a matrix L_z represents a finite difference scheme to calculate the derivative in z direction, i.e. ${\left({\mathbf{L}}_z\mathit{x}\right)}_j\approx \frac{\partial \mathit{\varphi }}{\partial z}{|}_{\mathit{r}={\mathit{r}}_j}$, and the diagonal matrix V represents the volumes of grid cells (V_ii being 1∕4 of the volume of all grid cells containing grid point r_i; see Sect. 2.5 for details), then

Other terms of the integral in Eq. (13) can be computed similarly. Using the notation of Eq. (14), with ${\left({\mathbf{L}}_{zz}\mathit{x}\right)}_j\approx \frac{{\partial }^{\mathrm{2}}\mathit{\varphi }}{\partial z^{\mathrm{2}}}{|}_{\mathit{r}={\mathit{r}}_j}$ and ${\mathbf{L}}_x,{\mathbf{L}}_y,{\mathbf{L}}_{xx},{\mathbf{L}}_{yy}$ defined similarly, we construct the precision matrix

with

To actually implement the calculation of Eq. (15) on an irregular grid one needs to have a triangulation for that grid, an interpolation algorithm (both described in Sect. 2.3). A method for calculating the spatial derivatives of the retrieved quantity at each grid point, needed to obtain ${\mathbf{L}}_i,{\mathbf{L}}_{ii},i=x,y,z$, is described in Sect. 2.4. The estimation of the physical parameters in Eq. (13) is discussed in Sect. 4.3. The errors introduced while constructing the precision matrix $S_{\mathrm{a}}^{-\mathrm{1}}$ are estimated in Appendix A. Applicability of methods presented in this section to 1-D and 2-D retrievals is briefly introduced in Appendix B.

2.3 Delaunay triangulation and interpolation

We aim to perform a retrieval on an arbitrary, finite grid. This section explains how vertically stretched Delaunay triangulation with linear interpolation is employed for that purpose.

To maintain generality and compatibility with arbitrary grids, we partition the retrieval volume into Euclidean simplexes (tetrahedrons in our 3-D case) with vertices at grid points. Many such partitions exist, but it is beneficial to ensure that the number of very elongated tetrahedrons, which increase interpolation and volume integration errors, are kept to a minimum. The standard technique used to achieve this is a Delaunay triangulation (Delaunay, 1934; Boissonnat and Yvinec, 1998), which maximises the minimum solid angle of any tetrahedron employed. Recall, also, from Sect. 2.2 that atmospheric quantities in UTLS and the stratosphere above tend to have much more variation vertically than they do horizontally within similar length scales. This difference can be quantified by comparing horizontal and vertical correlation lengths L_h and L_v. Hence, stretching the space vertically by means of coordinate transformation $\left(x,y,z\right)\mapsto \left(x,y,\mathit{\eta }z\right)$ with $\mathit{\eta }\sim L_{\mathrm{h}}/L_{\mathrm{v}}$ before constructing a Delaunay triangulation ensures that, on average, the amount of variation in the retrieved quantity in vertical and horizontal directions is similar. This, in turn, allows us to make use of the aforementioned benefits of Delaunay partitioning. The Computational Geometry Algorithms Library (CGAL, https://www.cgal.org, last access: 19 December 2018) was used to construct Delaunay triangulations for our grids.

In our implementation, we retrieve several quantities on the same grid, and L_h and L_v may differ for each of them. The precision matrix of Eq. (15) is evaluated separately for each quantity, so regularisation can still be calculated correctly as long as spatial structures of every quantity are adequately sampled by the grid. In practice, this means that as long as L_h∕L_v are of the same order of magnitude for all quantities, they can all be retrieved on one grid. Fortunately this is usually the case (see Sect. 4.3). For test retrievals in this paper, $L_{\mathrm{h}}/L_{\mathrm{v}}=\mathrm{200}$ for trace gas concentrations and $L_{\mathrm{h}}/L_{\mathrm{v}}=\mathrm{67}$ for temperature; therefore we set η=100. Retrieval results do not indicate systematic oversampling or undersampling in either horizontal or vertical directions and are not sensitive to a change in η.

In order to keep the computational costs low, we use a simple linear interpolation scheme: in each Delaunay cell we express an atmospheric quantity f at any point x as

If r_i, $\mathrm{0}\le i\le \mathrm{4}$ are the four vertices of the Delaunay cell, the (unique) constant gradient k can be obtained from a system of three linear equations:

This method ensures that the interpolated quantity is continuous and consistent with the volume integration scheme described in Sect. 2.5. Implementation is simple and fast, as any point inside a given Delaunay cell can be interpolated with data about that cell only. The gradient ∇f of the atmospheric quantity is, however, generally discontinuous at cell boundaries, making the interpolation unsuitable for direct use in spatial derivative evaluation.

2.4 Derivatives

In order to evaluate the cost function based on 3-D exponential covariance (Eqs. 13, 15), we need to estimate ∇ϕ and Δϕ at every grid point (as before, $\mathit{\varphi }=f-f_{\mathrm{apr}}$ is the departure of atmospheric quantity f from the a priori), i.e. construct the matrices ${\mathbf{L}}_i,{\mathbf{L}}_{ii},i=x,y,z$ of Eq. (15). This section describes the algorithm we use to achieve that on irregular Delaunay grids.

We begin by establishing some requirements that our derivative-estimation algorithm will have to meet. Let us write ${\mathbf{S}}_{\mathrm{a}}^{-\mathrm{1}}=\mathbf{A}+\mathbf{B}$, where A is a diagonal matrix representing the first integral of Eq. (13) (or first term of Eq. 10), and B represents the remaining terms. If B were an exact representation, it would be positive definite by construction, but this may not hold with ∇ϕ and Δϕ obtained numerically for a finite grid. If indeed an atmospheric state vector q exists so that q^TBq=0, then for every atmospheric state x

Hence, one of the states x±q would be favoured (have a smaller cost function) over x by the inverse modelling algorithm. Therefore, a scaled version of q would appear on any retrieval as noise. This particular type of noise would only be suppressed by the cost function term from ${\left(\mathit{x}-{\mathit{x}}_{\mathrm{apr}}\right)}^T\mathbf{A}\left(\mathit{x}-{\mathit{x}}_{\mathrm{apr}}\right)$. This guarantees only that x is not very far from the a priori, but does nothing to ensure that it is even continuous. This often results in retrievals with unphysical periodic structures (noise) and is unacceptable.

We deal with this problem by explicitly ensuring that B is positive definite. Consider the estimation of derivatives at a single grid point a∈ℝ³. ∇ϕ(a) and Δϕ(a) are numerically estimated from the values of ϕ in some (say m) grid points near a. We can write this as a map:

Now q^TBq=0, q≠0 implies that at every grid point a∈ℝ³ we have ∇ϕ(a)=0, Δϕ(a)=0 for the ϕ(a) corresponding to state q. Hence, such q will not exist if we require that D_a has trivial null space for all grid points a (i.e. $D_{\mathbf{a}}\left(\mathit{\varphi }\right)=\left\{\mathrm{0},\mathrm{\dots },\mathrm{0}\right\}$ only if ϕ(r_i)=ϕ(a), $\mathrm{1}\le i\le m$). The converse is not always true, so the trivial null space requirement is a stronger condition than absolutely necessary. It is, however, advantageous because it restricts the derivative-estimation algorithm at each grid point independently and hence can be implemented without computationally expensive operations on ${\mathbf{S}}_{\mathrm{a}}^{-\mathrm{1}}$ as a whole.

Since Eq. (20) will be used directly to construct ${\mathbf{S}}_{\mathrm{a}}^{-\mathrm{1}}$, which does not depend on atmospheric state x and hence on ϕ(a), D_a must be a linear map, so null space is a vector space of dimension $max\left\{\mathrm{0},m-\mathrm{6}\right\}$. Hence we need m≤6. We also want to avoid the derivative estimation to be underdetermined (m≥6), so we choose m=6; i.e. we aim to construct a linear map (Eq. 20) that estimates derivatives at grid point a using atmospheric quantity values at six points around a. This prevents us from using an interpolation for derivative calculations, since each interpolated value depends on atmospheric quantity values at multiple (in our case 4) grid points, and consequently, even if the same grid point is reused for several interpolated values, the total number of grid points employed exceeds six in any interpolation-based scheme we could come up with.

A solution satisfying the above criterion and able to deal with grid points with neighbours at irregular positions is based on polynomial fitting. For six grid points around a grid point $\mathit{a}=\left(x_a,y_a,z_a\right)$ denoted as ${\mathit{r}}_i=\left(x_a+x_i,x_a+y_i,x_a+z_i\right),\mathrm{1}\le i\le \mathrm{6}$, we write

We solve this as a linear system for the unknowns ${\mathit{\varphi }}_x,{\mathit{\varphi }}_y,{\mathit{\varphi }}_z,{\mathit{\varphi }}_{xx},{\mathit{\varphi }}_{yy},{\mathit{\varphi }}_{zz}$, which can then be directly identified with the derivatives we seek. Note that, as we required, the system can be solved by inverting a 6-by-6 matrix that only depends on r_i and not ϕ(r_i), so that we could construct the precision matrix (Eq. 15) to be used for every atmospheric state on a given grid.

The only remaining issue is selecting a suitable set of grid points r_i, so that the resulting matrix that needs to be inverted would not be singular or have a very high condition number. We found that this is much easier to achieve by imposing simple geometric criteria for point selection rather than algebraic conditions.

Let a be a point on Delaunay grid. The main idea of the algorithm we propose is to pick, for each spatial dimension d, a pair of points r₁, r₂ suitable for estimating the derivative in that direction. Points of such a pair should be sufficiently separated in d direction (i.e. projections of the vectors r_i−a in the d direction are sufficiently different), be on the opposite sides of a (unless a is at the boundary of the grid and this is not possible), and be reasonably close to a. We obtain the required six points by choosing one pair of points for each of the three dimensions. This can be implemented as follows.

Let $\mathit{a}=\left(x_a,y_a,z_a\right)$, and let us write the coordinates of other points on this grid as ${\mathit{r}}_i=\left(x_a+x_i,y_a+y_i,z_a+z_i\right)$ with $r_i=∥{\mathit{r}}_i-\mathit{a}∥$. Let $\mathrm{0}<\mathit{\beta }<\mathrm{1}$ and γ>1 be constants; suitable numeric values will be discussed further on.

1. For each Delaunay grid neighbour r_i of a, compute ${\mathit{\alpha }}_{ix}=x_i/r_i$, ${\mathit{\alpha }}_{iy}=y_i/r_i$, ${\mathit{\alpha }}_{iz}=z_i/r_i$.

2. If i, j such that $\left|{\mathit{\alpha }}_{ix}\right|>\mathit{\beta }$, $\left|{\mathit{\alpha }}_{jx}\right|>\mathit{\beta }$, α_ixα_jx<0 exist, then pick r_i, r_j for the derivative calculation. Else if k, l such that $\left|{\mathit{\alpha }}_{kx}\right|>\left|{\mathit{\alpha }}_{lx}\right|>\mathit{\beta }$, ${\mathit{\alpha }}_{kx}/{\mathit{\alpha }}_{lx}>\mathit{\gamma }$ exist, then pick r_l, r_k for derivative calculation.

3. Repeat step 2 for y and z dimensions. Each point can be only picked for one of the dimensions.

4. If neither condition of step 2 is satisfied for at least 1 dimension, repeat steps 1–3 with not only direct Delaunay grid neighbours of a, but also second neighbours (i.e. neighbours of neighbours). If this fails as well, do not calculate derivatives at the point a and simply assume them to be zero.

We found that the polynomial interpolation method is robust and a rather low value of β=0.3 was sufficient to produce suitable sets of neighbouring points (i.e. almost all 6-by-6 matrices resulting from such choice of points can be reliably inverted for polynomial fitting). In practice this means that the algorithm can find a suitable pair of points for dimension d unless all the Delaunay neighbours of the point a are located very close to the other two axes with respect to the point a.

2.5 Volume integration

Computing Eq. (13) requires a way to estimate volume integrals on Delaunay grid, which is described in this section.

Let a Delaunay cell have vertices r_i and face areas S_i for $\mathrm{0}\le i\le \mathrm{3}$ and volume V. Then, integrating an atmospheric quantity f as interpolated in Eq.  (17) over V we get

One can prove this by writing Eq. (17) with ${\mathit{r}}_i,\mathrm{0}\le i\le \mathrm{3}$ instead of r₀ and adding up the resulting four equations to obtain

Now choose Cartesian coordinates $\left(x^{\prime },y^{\prime },z^{\prime }\right)$, where ${\mathit{r}}_{\mathbf{0}}=\left(H,\mathrm{0},\mathrm{0}\right)$ and ${\mathit{r}}_i,\mathrm{1}\le i\le \mathrm{3}$ lie in the plane $x^{\prime }=\mathrm{0}$ and consider the x^′ component of the integral $\mathit{I}={\int }_V\left(\mathit{r}-\mathrm{1}/\mathrm{4}{\sum }_{i=\mathrm{0}}^{\mathrm{3}}{\mathit{r}}_i\right)\mathrm{d}V$:

As we can choose such coordinates for any vertex instead of r₀, the integral I has zero component in three linearly independent directions, so I=0. Then integrating Eq. (23) over V yields Eq. (22) and completes the proof.

We can now see that Eq. (22) gives an exact value to the volume integral of f, provided that f behaves exactly as prescribed by our interpolation scheme. Since the interpolation is linear, we may expect the values of volume integrals to be correct to the first order in cell dimensions. Since, generally, all data fed into the inverse modelling problem are interpolated beforehand, there is little reason to expect that a more elaborate and accurate volume integration scheme would improve the accuracy of retrievals on its own. A method based on Voronoi cells was also attempted and did not produce meaningfully different results.

4.1 The GLORIA instrument and data processing

The Gimballed Limb Observer of Air Radiance in the Atmosphere (GLORIA) is an aircraft-based infrared limb imager. It is a Michelson interferometer with a 2-D detector array, spectral coverage from 770 to 1400 cm⁻¹ and spectral sampling of up to 0.0625 cm⁻¹ (Riese et al., 2014; Friedl-Vallon et al., 2014). It is carried aboard the Russian M55 Geophysica or German HALO research aircraft and is generally intended to look to the side of the aircraft, but can pivot horizontally so that the angle between aircraft heading and observation direction can be changed from 45 to 135^∘. The limb observation set-up naturally allows for high vertical resolution (up to approximately 200 m) and, with a proper trade-off between spectral resolution and measurement time, a high density of vertical profiles sampled along-track. The pivoting allows us to look at the same air mass from several points along the flight path and use tomography to improve resolution in the horizontal direction perpendicular to the flight path. The best results with tomographic retrievals are achieved, however, when the aircraft flies in a path close to circular (hexagonal flight paths are used in practice, so that the aircraft could fly straight most of the time) of around 400 km in diameter, and GLORIA observes the air masses in the middle from many directions (Ungermann et al., 2011). A horizontal resolution down to 25 km in both directions can be achieved in this way. As with any limb observer, GLORIA can only provide tangent point data about air masses at and below aircraft flight altitude, which is limited to approximately 20 km for the M55 Geophysica and 15 km for HALO. For more detail on GLORIA, refer to Friedl-Vallon et al. (2014), for example.

The implementations of the algorithms described in the Sects. 2 and 3 were integrated into the Jülich Rapid Spectral Simulation Code version 2 (JURASSIC2). The forward model of atmospheric radiance used in this code employs a forward model for atmospheric radiation based on radiances obtained from the emissivity growth approximation method (Weinreb and Neuendorffer, 1973; Gordley and Russell, 1981) and the Curtis–Godson approximation (Curtis, 1952; Godson, 1953). For more information about JURASSIC2, refer to Hoffmann et al. (2008); Ungermann et al. (2011); Ungermann (2013). The performance of the new algorithms, both in terms of output quality and computational cost, is evaluated in the rest of this section.

Figure 1Panel (a) shows the ECMWF temperature data used to generate the synthetic measurements (true temperature). Panel (b) shows the temperature residual: the difference between the true temperature and the smoothed temperature field used as an a priori. In the first column, black lines represent the contours of a priori temperature values. In columns 2 and 3, the black line is the flight path, the thick grey line is the position of the cut shown in the first column and the thin grey line is the topography (Icelandic coastline visible).

Download

4.2 Retrieval set-up and test data

A hexagonal flight pattern intended for tomography was realised on 25 January 2016, as part of flight 10 of the POLSTRACC measurement campaign. The flight path of the HALO research aircraft contained a regular hexagon with a diameter (distance between the opposite vertices) of around 500 km over Iceland, which allowed for a high spatial resolution retrieval in the central part of the aforementioned hexagon (Fig. 1). The flight altitude throughout the hexagonal flight segment remained close to 14 km. For detailed information about this flight refer to Krisch et al. (2017).

The test case used in this paper is based on the actual aircraft path, measurement locations, spectral lines used for retrieval and meteorological situation during this flight. Synthetic measurement data were used instead of real GLORIA observations: the forward model of JURASSIC2 was employed to simulate the observed radiances in the atmospheric state given by the ECWMF temperature and WACCM model data for trace gases. Simulated instrument noise was subsequently added to those radiances. This set-up allows us to use the model data as a reference (the “true” atmospheric state) for evaluation of retrieval quality. The same set of simulated measurements was used for all test retrievals described in this paper. These measurements were obtained by running the forward model on a very dense grid (about twice as dense in each dimension as those of the densest test retrievals). This was done to ensure that the discretisation errors in the simulated measurements would be minimal and would not favour any retrieval (as could happen if they were generated on the same grid). An evaluation of forward model errors on different grids is presented in Appendix C.

The retrieval derives temperature and volume mixing ratios of O₃, HNO₃, CCl₄ and ClONO₂. Temperature and CCl₄ concentration are derived in the altitude range from 3 to 20 km; O₃, HNO₃ and ClONO₂ are retrieved from 3 to 64 km. The latter three trace gases have larger volume mixing ratios above flight level and hence significant contributions to measured radiances from these high altitudes. Interpolated WACCM data were used as an a priori for all trace gases. A priori data for temperature were obtained by applying polynomial smoothing to ECMWF data. The retrieval is hence not supplied with the full temperature structure that was used to simulate the measurement. Therefore, the agreement between the retrieval result and the “true” atmospheric state, upon which the simulated measurements are based, cannot be achieved by overregularisation.

Figure 2Panels (a) and (b): horizontal cuts of retrieval grids. Dots show grid points; blue lines indicate Delaunay cell boundaries. Panel (a) shows the full initial grid for retrievals A, B and C (all horizontal layers of this grid are identical), (b) shows the thinned-out grid for retrieval D (cut at 11 km). Panel (c) is the grid for retrieval D. Radial distance is defined as distance to the vertical line at the centre of the hexagonal flight path (66^∘ N, 15^∘ W). Points within the inner core (1000 km ×1000 km ×20 km) are shown in blue, other points are in red, and the shaded rectangle shows the position of the horizontal cut in panel (b).

Download

4.3 Correlation length estimation

Table 1Retrieval parameters.

Regularisation weights α₀, α_h, (as in Eq. 2) are used for retrieval A, the standard deviation σ (see Eq. 15) used for retrievals B-D.

Download Print Version | Download XLSX

The standard deviation σ and correlation of the lengths L_h and L_v in Eq. (13) can be obtained from statistical analysis of in situ measurement or model data. Rigorous derivation of these parameters for each retrieved trace gas and temperature would require detailed analysis of in situ measurement databases and is, therefore, out of scope of this paper. However, the following rough approximation proved sufficient to improve the retrieval results compared to the old regularisation scheme.

3-D tomographic retrievals are most useful for those trace gases that have high vertical gradients and complex spatial structure. The spatial correlation lengths of concentrations of such gases are mostly determined by stirring, mixing and other dynamical processes and are therefore similar. Furthermore, Haynes and Anglade (1997) and Charney (1971) estimated the dynamically determined aspect ratio to be α≈200–250 in the lower stratosphere, and we expect less than that in the upper troposphere. Hence, we will assume $L_{\mathrm{h}}/L_{\mathrm{v}}=\mathrm{200}$ for all trace gases.

An analysis of spatial variability of some airborne in situ measurements was performed by Sparling et al. (2006). For ozone concentration χ(s ,z), where s is the horizontal location and z is the altitude, they define fractional difference parameters:

and provide the dependence of the standard deviation σ(Δ_r) on the horizontal separation r. If the typical spatial variation in ozone concentration is significantly smaller than the mean ozone concentration 〈χ〉, one can approximate ${\mathrm{\Delta }}_r\approx \left|\mathit{\chi }\left(s+r\right)-\mathit{\chi }\left(s\right)\right|/\langle \mathit{\chi }\rangle$. Then, using the covariance relation (8) and various properties of the normal distribution, we get ${\mathit{\sigma }}^{\mathrm{2}}\left({\mathrm{\Delta }}_r\right)=\mathrm{2}{\left(\mathit{\sigma }/\langle \mathit{\chi }\rangle \right)}^{\mathrm{2}}\left(\mathrm{1}-\mathrm{2}/\mathit{\pi }\right)\left(\mathrm{1}-\exp \left(-r/L_{\mathrm{h}}\right)\right)$. By comparing this to the experimental results of Sparling et al. (2006), we estimate L_h=200 km. Using the aspect ratio argument above, we also set L_v=1 km.

The spatial structure of retrieved temperature differs from that of the trace gases. It is determined not only by mixing and stirring, but also radiative processes and gravity waves. Radiative transfer tends to erase some of the fine vertical structure determined by isentropic transport; hence one should use a greater vertical correlation length than in the case of trace gases. Gravity waves have a variety of length scales, but due to the finite resolution of the instrument not all of them can be retrieved. Following the gravity wave observational filter study for the GLORIA instrument in Krisch et al. (2018) we can estimate the lowest observable gravity wave horizontal wavelength to be λ_h∼200 km and the lowest vertical wavelength to be λ_v∼3 km. λ_h coincides with the dynamically determined correlation length L_h=200 km, so this value should be appropriate for temperature as well. The λ_v limit is not only determined by instrument resolution, but is also related to the fundamental properties of gravity waves. Preusse et al. (2008) gives the following expression for the gravity wave saturation amplitude

where $\overline{T}$ is the mean temperature and N is the Brunt–Väisälä frequency. In typical lower stratosphere conditions and when λ_v=3 km, Eq. (33) gives ${\widehat{T}}_{\max }\approx \mathrm{4}\mathrm{K}$. Since gravity waves are usually not saturated and considering typical measurement errors (see Sect. 4.6), waves of significantly shorter wavelengths would probably not be observed. We hence set L_v=3 km. Time series of measurements at a fixed point in the atmosphere are then used to approximate the standard deviations σ for all retrieved quantities (Table 1). Our choice of σ value for temperature, in particular, may seem rather low. This is related to the choice of a priori: smoothed ECMWF data were used instead of climatologies that would have been more typical in this case. Such an a priori can be expected to match the large-scale structures of real temperature more closely, thus reducing the forward model errors in poorly resolved areas and improving the results. It is particularly useful for resolving fine structures, such as gravity waves, which were the main scientific interest of the measurement flight used as a basis for test retrievals (Krisch et al., 2017). A likely closer match between the a priori and retrieval thus requires a lower σ value, that, in this set-up, represents expected strength of gravity wave disturbances, rather than the full thermal variability of the atmospheric region in question.

Table 2Grid densities in different regions for retrieval D in kilometres.

Download Print Version | Download XLSX

4.4 Test retrievals

The test data described in the previous section were processed with several different regularisation set-ups and retrieval grids. Here we present the results for four different runs.

Firstly, as a reference, the latest version of JURRASIC2 without the implementation of any new algorithms described in this paper was used (retrieval A). It uses a rectangular grid that covers altitudes from 1 to 64 km. The vertical spacing of the grid is 1 km between the altitudes of 1 and 5 km, 125 m between altitudes of 5 and 20, 1 km between altitudes of 20 and 25 and 4 between 28 and 64 km. As required for the rectilinear grid, the number and distribution of points is the same in each altitude (Fig. 2a). For this reference run, a first-order regularisation (Eq. 1) with trilinear interpolation (Bai and Wang, 2010) was used. The regularisation weights, as defined in Eq. (2), are given in Table 1. These are typically tuned ad-hoc, i.e. adjusted by trial and error until optimal retrieval results can be achieved, starting from the default values obtained as in Hansen and O'Leary (1993). While tuning, the retrieval results are validated against the atmospheric state used to generate the synthetic measurements. When using real observations, the results are validated against in situ measurements and model data.

Figure 3Panels show the difference between retrieved temperature and true temperature of the simulated atmosphere. Rows A–D show the results of the respective retrievals. In the first column, black lines represent the contours of the a priori temperature value. In columns 2 and 3, the black line is the flight path, the thick grey line is the position of the cut shown in the first column and the thin grey line is the topography (Icelandic coastline visible).

Download

Then, to evaluate the performance of second-order regularisation, retrieval B was performed. The grid and interpolation methods were identical to retrieval A, but the regularisation was replaced with the second-order scheme from Eqs. (13) and (15) and correlation lengths derived in Sect. 4.3 from in situ observations. No subsequent tuning of these parameters was performed.

We compare the quality of different retrievals by inspecting the differences between the retrieved temperature and temperature used to generate the synthetic measurements (“true” temperature). These differences are shown on horizontal and vertical slices of the observed atmospheric volume in Fig. 3. Figure 4 shows the variation of the retrieved and true temperatures on horizontal lines in order to visualise the response of retrievals to different spatial structures in true temperature.

Figure 4Retrieved temperature at different altitudes in the vertical cut shown in Figs. 1 and 3. Line T represents the true temperature of the simulated atmosphere, A–D are the respective retrievals.

Download

A comparison of Fig. 3A and B shows the effect of the new regularisation (Eq. 13). The new algorithm (retrieval B) demonstrates better agreement with the true temperature. It shows markedly less retrieval noise in the central area, as well as better results (less effect of a priori) just outside it, at the edges of the vertical cut shown in Fig. 3. Considering also that, unlike A, the regularisation strength of B was not manually tuned for this test case in particular, we can conclude that the new regularisation performed better in this test case.

The third retrieval (C) was performed on the same set of grid points as B and with the same input parameters, but the grid was treated as irregular; i.e. a Delaunay triangulation was found. The algorithms described in Sect. 2 were used for regularisation. A horizontal cut of the Delaunay triangulation for this grid is shown in Fig. 2a. By comparing B and C one can evaluate the agreement between the old interpolation, derivative calculation and volume integration techniques for rectilinear grids and their newly developed irregular-grid-capable alternatives. The temperature fields from retrievals B and C are very similar (Fig. 3), which confirms that the new Delaunay-triangulation-based code gives consistent results when run on a rectilinear grid.

Finally, retrieval D was performed with the same input parameters and methods as C, but using only the subset (22 %) of the original grid points to reduce computational cost. These points were chosen so that grid density would be highest in the volumes best resolved by GLORIA and sparser where little measurement data are available. In particular, a set of axially symmetric regions around the centre of original retrieval (i.e. symmetric with respect to vertical axis at 66^∘ N, 15^∘ W) was defined as shown in Table 2. Points were removed from the original grid to achieve the specified vertical and horizontal resolutions for each region, resulting in the grid shown in Fig. 2. Results changed little compared to retrieval C. Most of the minor differences between C and D occur near the edges of the hexagonal flight path (Fig. 3). Even though the grid spacing for retrieval D is still the same as C in some these areas, larger grid cells nearby have some effect on their neighbours. Note also that C and D retrievals tend to differ in the areas where the agreement between each of them and the true temperature is relatively poor and diagnostics predict lower accuracy (Fig. 5), so D can be considered to be as good as C in well-resolved areas. Unphysical temperature structures far outside the measured area, like the one partially seen at 64^∘ N, 20^∘ W, are not unique to this retrieval. They just occur a little closer to the flight path in this case due to the coarser grid outside the hexagon.

Figure 5Panels B and D show the Monte Carlo estimate of the total temperature error the retrievals B and D would have if the measurements were real. Panel Ds shows an error estimate for the synthetic retrieval D (this assumes that unretrieved quantities, e.g. pressure, are known exactly). In columns 2 and 3, the black line is the flight path, the thick grey line is the position of the cut shown in the first column and the thin grey line is the topography.

Download

4.5 Computational costs

JURASSIC2 finds an atmospheric state minimising Eq. (1) iteratively. In each iteration it first computes a Jacobian matrix of the forward model and hence obtains a Jacobian of the cost function (1), then minimises the cost function using the Levenberg–Marquardt (Marquardt, 1963) algorithm and by employing conjugate gradients (CG) (Saad, 2003, Hanke, 1995) to solve the linear systems involved. The forward model is then run to simulate the measurements for the new iteration of the atmospheric state.

All calculations were performed on the Jülich Research on Exascale Cluster Architectures (JURECA) supercomputer operated by the Jülich Supercomputing Centre (Jülich Supercomputing Centre, 2016). Each retrieval was executed in parallel on four computational nodes, each with 24 CPU cores (twin Intel Xeon E5-2680 v3 CPUs) and 128 GB of system main memory.

Table 3 shows the time required for different parts of calculation averaged over three identical runs (although <1 % variation was observed). The “other calculations” field mostly includes the time required by input and output, but also preparatory steps such as the Delaunay triangulation in the case of retrievals C and D.

The results in Table 3 reveal the impact of the new methods on computation times. The new, more complex regularisation results in slower convergence in the conjugate gradient calculation, making retrieval B about 6 % more expensive than A. The Delaunay triangulation, however, decreases Jacobian calculation costs by roughly 25 %, even the one used on the same grid as the old approach. This has to do with the fact that interpolation based on Delaunay triangulation only requires four values of the interpolated quantity at neighbouring points, and linear interpolation in three directions requires eight, thus resulting in more non-zero terms in the Jacobian. The time required for the initialisation of the Delaunay interpolation and derivatives, including construction of the matrix ${\mathbf{S}}_{\mathrm{a}}^{-\mathrm{1}}$, increases rather fast with grid size with current implementation, contributing to “other calculations” of retrieval C. Further optimisation of this step is possible and will be implemented in the future. In total, the new retrieval (C) is 9 % faster then the old approach (A) while using the same grid.

Table 3Computational cost of retrievals (calculation times given in minutes).

Download Print Version | Download XLSX

Retrieval D gives a further 45 % reduction of the computation time by removing 82 % of points in the grid. A speed-up of such an order was expected. Only a minority of grid points contribute to a forward model (and hence Jacobian) calculations, and JURASSIC2 code is heavily optimised for sparse matrix operations. Most of the grid points contributing to the forward model are in the areas well resolved by the instrument and thus cannot be removed from the grid without compromising data products. Therefore, the actual savings in computation time mostly come from removing points above a flight level and choosing a more appropriate non-rectangular shape for the densest part of the grid when measurement geometry requires it. Some of the points above the flight path must be included since infrared radiation from there is measured by the instrument, but very little information about atmosphere far above the flight level can be retrieved, so very low grid densities are sufficient (Table 2).

4.6 Diagnostics

We follow Ungermann et al. (2011) and use Monte Carlo approach to obtain diagnostics, but replace Cholesky decomposition with the method developed in Sect. 3. Results for retrievals B and D are presented in Fig. 5B and D. They show the total error estimate for temperature based on the sensitivity of temperature value to (simulated) instrument noise and all other retrieved atmospheric quantities as well as the expected error of non-retrieved quantities (e.g. pressure). These results tell us what temperature errors we could expect if the retrievals in question were performed with real data. Panel Ds shows an error estimate based on (simulated) instrument noise and uncertainties of retrieved parameters only. In the case of a synthetic retrieval, unretrieved parameters, upon which the simulated measurements are based, are known exactly. Therefore, panel Ds is a direct Monte Carlo prediction of the temperature deviation seen in Fig. 3D, and the magnitude of this deviation is indeed similar.

These results confirm that temperature can only be reliably retrieved within the flight hexagon or very close to the actual path. Also, the error estimates are higher in the central area of hexagon, close to flight altitude, in agreement with test results (compare column 3 of Figs. 3 and 5).

Monte Carlo retrieval results provided here were verified by estimating the temperature error at one grid point using a more direct approach. Rodgers (2000) proves that the retrieved atmospheric state x can be written in terms of the a priori state x_a, the (unknown) true state x_t and measurement error ϵ as

where

Here, D denotes the Jacobian operator and all other notation as in Eq. (1). The so-called gain matrix G fully describes the effect of measurement uncertainties ϵ on the retrieved state x. Computing G would be prohibitively expensive due to the required matrix inversion, but note that the i'th element of x only depends on the i'th line of the matrix G, so one can use Krylov subspace methods (e.g. conjugate gradients) to solve Mr_i=e_i for the i'th row ${\mathit{r}}_i^T$ of the inverse matrix M⁻¹. The i'th row of G can then be obtained as ${\mathit{r}}_i^T{\left(\mathbf{DF}\left(\mathit{x}\right)\right)}^T{\mathbf{S}}_{\mathit{ϵ}}^{-\mathrm{1}}$.

This approach allows for error estimation for a small number of grid points in a reasonable amount of time. We chose one point at the centre of the hexagon (66^∘ N, 15^∘ W, 11.75 km altitude) for retrieval B. The total temperature error was found to be 1.27 K and agrees with the Monte Carlo result (1.25 K) at this point within 2 %, which is an excellent accuracy for an error estimate.