Characterization of model errors in the calculation of tangent heights for atmospheric infrared limb measurements

We review the main factors driving the calculation of the tangent height of spaceborne limb measurements: the ray-tracing method, the refractive index model and the assumed atmosphere. We find that commonly used ray tracing and refraction models are very accurate, at least in the midinfrared. The factor with largest effect in the tangent height calculation is the assumed atmosphere. Using a climatological model in place of the real atmosphere may cause tangent height errors up to ± 200 m. Depending on the adopted retrieval scheme, these errors may have a significant impact on the derived profiles.


Introduction
Inversion algorithms for atmospheric limb measurements from space usually retrieve profiles on pressure coordinates due to the significant uncertainty of the line of sight.Accurate knowledge of the line of sight is however needed to reconstruct the altitude grid of the retrieved profiles, as it may be necessary for comparison to correlative measurements (such as those obtained from lidar) that are intrinsically represented on an absolute altitude grid.Spectral measurements also contain information on the instrument viewing direction; however, the spectral resolution and the signalto-noise ratio are often insufficient to determine accurate estimates of the line of sight.The accurate calculation of the line of sight from the engineering estimates of the instrument pointing angles and instrument position relies, however, on the accuracy of both the ray-tracing algorithm and the model used for atmospheric refraction.In this work, we compare the accuracy of a few ray-tracing and atmospheric refraction models using the tangent height error as a quantifier.The tests presented are based on measurements of the Michelson Interferometer for Passive Atmospheric Sounding (MIPAS, Fischer et al., 2008) that successfully operated on board of the polar satellite ENVISAT in the time frame from April 2002 to April 2012.While based on the setup of a specific instrument, the algorithms investigated are applicable to the general problem of ray tracing in an inhomogeneous transparent atmosphere.

Ray tracing
The propagation path of electromagnetic rays through an inhomogeneous medium can be deduced from the eikonal equation (Born and Wolf, 1970): where x is the position vector, n(x) the refractive index and φ(x) is the so called eikonal function.This equation can be derived directly from the first-order Maxwell equations or from the second-order wave equations for either the electric or the magnetic field.The only simplifying assumption used in this derivation is that n(x) varies slowly with respect to the wavelength of the electromagnetic field.Of course, for the propagation of mid-infrared radiation in the Earth's atmosphere, this hypothesis is verified with very high accuracy.The surfaces φ(x) = constant are the geometrical wave fronts.Therefore, the ray direction is parallel to ∇φ(x).Let p(s) be the ray path, parametrized with the arc parameter s.
Published by Copernicus Publications on behalf of the European Geosciences Union.
We can write dp(s) ds = ∇φ(p(s)) |∇φ(p(s))| . (2) After some algebraic manipulation, using Eq. ( 1) we get the differential form of the light rays equation (Born and Wolf, 1970): This is a vectorial second-order differential equation that permits one to derive the full ray path across an inhomogeneous medium, if n(x) and the boundary conditions are known.From this equation, several ray-tracing methods can be derived with different tradeoffs between accuracy and computational speed.In this work we consider the following three ray-tracing methods: -direct numerical solution of Eq. ( 3), -tangential displacement method (Hase and Höpfner, 1999), -iterative Snell's law (Thayer, 1967;Hobiger et al., 2008).
In Thayer's implementation of Snell's law the atmosphere is assumed horizontally homogeneous.This implementation is one of the fastest ray-tracing methods, however, if the horizontal variability of the atmosphere is taken into account, this method is not adequate.For horizontally varying atmosphere, the level lines of n(x) do not coincide with the altitude levels.Thus, in the Thayer implementation, the calculation of the refraction angle with Snell's law is based on a wrong hypothesis.In our tests we found that in several synthetic but realistic atmospheric conditions, with strong horizontal gradients of pressure and temperature, the method produces a wrong ray path, partly following Earth's curvature.In Hobiger et al. (2008), a refined approach of Thayer's method is proposed, removing the horizontal homogeneity assumption at the expense of significantly increased computational complexity.
The tangential displacement method (henceforth referred to as TD) is an iterative approach for the solution of the eikonal equation, using an approximation to avoid the calculation of the second derivatives of p(s).
For the direct numerical solution of the eikonal equation we implemented a multi-step predictor-corrector method (henceforth referred to as EIK) using the two-step Adams-Bashforth formula for the predictor and the BDF2 formula for the corrector (Isaacson, 1994).The shape of the ray path, however, suggests the use of an adaptive step length based on the second derivatives of p(s) that are linked to the local ray curvature.These derivatives are easily obtained from the numerical solution of the eikonal equation.Thus, we also implemented this adaptive method (henceforth referred to as AEIK) while still maintaining the property that in each atmospheric layer the step is fixed.This is an efficient choice in view of the implementation of the Curtis-Godson (CG) integrals for the calculation of the radiative transfer in a horizontally varying atmosphere.In fact, since the curvature of the ray path and the atmospheric variability of the quantities to be integrated are relatively small, it is possible to a priori estimate the step size for the numerical calculation of the CG integrals within a prescribed accuracy.With this approach the same set of nodes can be used for the calculation of all CG integrals, and the iterative refinement of the integration step size can be avoided.
All the implemented methods can be applied to a threedimensional ray tracing.Our implementation is however planned for inclusion in the ESA retrieval model for the Michelson Interferometer for Passive Atmospheric Sounding (MIPAS) routine data processing (Ridolfi et al., 2000;Raspollini et al., 2006Raspollini et al., , 2013)), which assumes that all limb scans lie in the orbit plane.In a forthcoming evolution of this algorithm, the atmosphere will be represented with a 2-D discretization in the orbit plane and will be considered constant in the direction perpendicular to this plane.Under this assumption all the ray paths lie in this plane, therefore we implemented a 2-D ray-tracing scheme.
A good online comparison of refraction formulas can be found in Young (2011).
The Barrel-Sears empirical formula has been used for a long time for atmospheric infrared applications.We include this formula in our tests mainly for historical reasons.The version implemented in our code is where p is the total air pressure expressed in hPa, T is the temperature in Kelvin, λ the wavelength in µm and p w is the water vapor partial pressure in hPa.
The simplified Edlén formula is the model currently implemented in the ESA retrieval code for routine MIPAS data inversion.The formula implemented in our code is with the constants p 0 = 1013.25 hPa, T 0 = 288.16K and c 0 = 0.000272632.This formula clearly does not model the dependence of the refraction on the wavelength and the water vapor amount.
Ciddor's formula models the refractive index as a function of wavelength, pressure, temperature, water vapor and carbon dioxide content.The formula was originally tested with experimental data extending only up to 1.7 µm; however, the work of Mathar (2004) suggests that its accuracy is on the order of 10 −6 also up to 25 µm, i.e., over the whole spectral region covered by MIPAS observations (4.1-14.9µm).We implemented Ciddor's formula following the original paper (Ciddor, 1996).

Atmospheric models
Whatever refraction model is chosen, the refractivity depends on pressure, temperature and, possibly, water vapor and carbon dioxide.Thus, the ray tracing depends on the assumed atmosphere.To evaluate the impact of the selected atmosphere we considered the following models: -the US Standard Atmosphere, 1976(US Gov., 1976), -the IG2 atmosphere (Remedios et al., 2007), -the atmosphere retrieved in a previous processing version of MIPAS data (Raspollini et al., 2013), -atmospheric refractivity profiles determined from co-located Radio Occultation (RO) measurements (Schwärz et al., 2012).
The US Standard Atmosphere (US76), together with the simplified Edlén formula for refraction, is the model currently adopted by ESA in MIPAS level 1b data processing (Kleinert et al., 2007) for determining the tangent height of the limb measurements, starting from the instrument position and pointing angles.The US76 atmosphere does not include any model for the horizontal or seasonal variabilities.
The IG2 climatological database is a collection of atmospheric profiles used as initial guess (IG) or assumed profiles in MIPAS routine level 2 retrievals (Remedios et al., 2007).The IG2 database includes vertical profiles of pressure, temperature and volume mixing ratio (VMR) of constituents relevant for MIPAS data processing.The profiles are tabulated as a function of year, season and latitude band.In our tests, we also linearly interpolate the profiles in latitude, in between the tabulated values.This interpolation makes the atmosphere smoothly change with latitude, avoiding abrupt changes at the edges of the latitude bands.
The last two atmospheres rely on experimental data.The tests with RO refractivity measurements considered in this work are limited to the MIPAS orbit 43442 acquired on 21 June 2010.In this orbit, there are 16 MIPAS limb scans for which a co-located RO measurement exists within 300 km and 3 h.

Results
The calculated ray path depends on the ray-tracing method, the refractive index model and the assumed atmosphere.To evaluate the impact of each of these factors, we use the error on the calculated tangent height of the limb measurements as a quantifier.The MIPAS level 1b data files provide the geolocation of the tangent points of the limb measurements.These are determined using the position and attitude of the satellite via a ray-tracing algorithm that uses the US76 atmosphere and the simplified Edlén formula for refractivity.Since we have no access to the algorithm details, however, we are not able to reproduce exactly the tangent height values reported in the level 1b data files.To maintain full consistency with the level 1b data, instead of starting the ray tracing from the satellite position and attitude, we use the following approach.We start from the level 1b tangent point geolocation and use the same assumptions of the level 1b processor to back-calculate the latitude and the slope of the ray path at the intersection with the atmospheric boundary, fixed at 120 km (outgoing path).Then we reverse the ray tracing (incoming path) and recalculate the tangent point.The assumptions used for the incoming path can be different from those used for the outgoing path.The difference in height between the original and the recalculated tangent point characterizes the impact of the different assumptions used for the incoming path.
All the tests reported below were carried out using the MI-PAS measurements acquired on 12 orbits (3 for each season) in the year 2010.We find that the estimated errors are correlated with latitude and season; however, their maximal amplitude changes very little in our sample of orbits.For this reason, here we only present the results of orbit 43442, for which correlative RO refractivity measurements are also available for some scans.
In order to assess the accuracy and computational efficiency of the ray-tracing methods described in Sect.2, we use the following test.For each limb scanning, we trace the incoming path with the same hypotheses used to calculate the outgoing path, i.e., US76 atmosphere and simplified Edlén refractivity model.Thus, the difference in the recalculated tangent heights is entirely due to the numerical error of the ray-tracing methods.Of course, the accuracy of each method depends on the adopted step length.The smaller the step length, the more accurate the solution; however, smaller step lengths require more computing time.To quantify the accuracy of a method, we use the summation on the limb scans of the whole orbit of the absolute values of the differences in the recalculated tangent heights.
In Fig. 1 we show the tradeoff between accuracy and computing time for the different methods described for various step lengths.In the case of the AEIK method, the values reported in the figure are the initial step lengths, which are then adapted by the method itself.For sufficiently small step lengths, all the three considered methods are very accurate.The AEIK method is the most efficient; it however requires With the fixed US76 atmosphere, we then tested the impact of the refractive index model used for the incoming path.In the case of a horizontally homogeneous atmosphere such as the US76, the impact of changing the refraction model can be estimated by the change in tangent height calculated from the constant value of the product r n(r) (where r is the distance of the tangent point from the Earth's center).However, we preferred to determine the tangent height error via the full ray-tracing algorithm that we had already implemented.Differences between original and recalculated tangent heights are always much smaller than the accuracy required for tangent heights (50-100 m).The simplified Edlén and Ciddor formulas turn out to be in very good agreement, providing differences in tangent heights of less than 21 cm.The old Barrel-Sears formula provides slightly more differing results compared to the other two refraction models (within 2.3 m).The same test repeated with the retrieved atmosphere instead of the US76 atmosphere does not change the magnitude of the differences found.In our subsequent tests, we use the Ciddor formula that is considered the most accurate (Young, 2011), and we use Edlén's formula as a backup option, should the water profile be unavailable.
Finally, with the same strategy we studied the impact of the assumed atmosphere on the ray tracing.As expected, we found that this is the assumption with the largest impact on the calculation of the height of the tangent points.In Fig. 2, we show the differences (color scale) in meters between the original and recalculated tangent heights, as a function of the orbital coordinate (horizontal axis) and height (vertical axis).We only show tangent heights less than 20 km, since for higher tangent heights the difference is less than 10 m.In Fig. 2a, the assumed atmosphere for the incoming path is the IG2 atmosphere, in Fig. 2b the retrieved atmosphere.We note that above the troposphere the two panels show good agreement.In the troposphere, where the atmospheric variability is larger, there are some differences.The larger differences are observed for the retrieved atmosphere and the lowest tangent heights where values up to 200 m are achieved.
For the 16 limb scans of orbit 43442 for which a co-located RO measurement is available, we repeated the test of Fig. 2 using the RO refractivity profiles for the calculation of the incoming path.For this test, the high vertical resolution RO refractivity profiles were first interpolated to the common vertical grid of the profiles used in the ray tracing (1 km).Moreover, in this case we assumed a horizontally homogeneous atmosphere.The results of this test are very similar to those obtained with the retrieved atmosphere (see Fig. 2b); therefore, we do not show the related map.
Note that the error in the tangent height shown in Fig. 2b includes the contributions due to the vertical variation of the atmosphere and to the horizontal gradients.In order to quantify these contributions, we repeated the test of Fig. 2b with the retrieved atmosphere but assumed the horizontal homogeneity in the ray tracing of each scan.We found that, on average, the two contributions are correlated, 85 % of the error is due to the vertical structure of the atmosphere and the remaining 15 % is due to the horizontal gradients.
In order to assess the influence of the newly calculated tangent heights on the level 2 products, we retrieved orbit 43442 starting from the recalculated tangent heights using the ESA operational algorithm (Ridolfi et al., 2000;Raspollini et al., 2006Raspollini et al., , 2013)).On average, the differences in the retrieved profiles are much smaller than the noise error, with no real improvement in the final χ 2 of the fit.However, locally, in the scans with large height corrections, the differences can be on the order of the noise error.
The small size of the differences in the retrieved profiles is due to the ability of the ESA retrieval code to adjust the pressure at the tangent points and recalculate tangent height increments using the hydrostatic equilibrium.Like the ESA operational code, several MIPAS inversion algorithms (Dudhia et al., 2005;Stiller, 2000) retrieve tangent height or pressure to compensate for the ray-tracing errors.On the other hand, due to the already high code complexity, the 2-D geofit multi-target retrieval (GMTR) retrieval algorithm (Carlotti et al., 2006) assumes correct instrument pointing.To show the compensation effect mentioned above, in the ESA code we forced the retrieval program to keep constant the tangent pressure during the iterations.In this case, the retrieved profiles are more sensitive to the used tangent height values.In Fig. 3, we show the influence on the temperature profile, retrieved using the original (T L1b ) and recalculated (T RET ) tangent heights.Figure 3a refers to the standard ESA retrieval algorithm, while in Fig. 3b we impose a constant tangent pressure.In both plots the average profile of T RET − T L1b is shown (blue curve), together with its standard deviation (gold curve) and the average noise error of the T RET (magenta curve) and T L1b (red curve) profiles.The magenta and red curves in both panels are almost identical.Note that both the average difference and the standard deviation are larger in Fig. 3b.At the lowest altitudes there are significant values of T RET − T L1b , with differences as large as 5 K, correlated with the largest tangent height corrections plotted in Fig. 2b.Since the tangent pressure is kept fixed in the retrieval of Fig. 3b, the algorithm can not use this parameter to adjust the tangent heights.Therefore, temperature is used, via hydrostatic equilibrium, to compensate for the error in tangent heights.

Conclusions
We analyzed the main factors driving the calculation of the tangent heights of spaceborne limb measurements.We found that the factor with largest effect in the tangent height calculation is the assumed atmosphere.Using a climatological model in place of the real atmosphere may cause tangent height errors up to ±200 m.This error is on the same order of magnitude of the most recent estimates of the error on MIPAS tangent height due to uncertainties on the pointing angles and satellite attitude.In MIPAS retrievals, this inaccuracy causes a temperature error on the order of the noise error if the tangent height is adjusted by fitting the tangent pressure (this is the case of the ESA algorithm).However, if the retrieval assumes a fixed tangent pressure, the inaccuracies in tangent heights may cause temperature differences locally exceeding 4-5 K.

Figure 1 .
Figure 1.Efficiency of the tested ray-tracing methods, for step sizes of 0.1, 0.25, 0.5 and 1.0 km.

Figure 3 .
Figure3.Average difference between temperature retrieved with recalculated (T RET ) and original (T L1b ) tangent heights, and standard deviation of the differences.Errors due to measurement noise on T RET and T L1b (the magenta profiles are obscured by the overlapping red profiles).Standard (a) and constant tangent pressures (b) retrieval cases.