2.4.4.1.3 The retrieval modules
2.4.4.1.3.1 General features of the
adopted approach
In contrast with the data processing
of already flown operational limb
sounding instruments (e.g. on the
Nimbus7 and UARS Satellites),
which have been mainly of
radiometric type, MIPAS data
analysis will be confronted with the
exploitation of broad band and high
resolution spectral measurements
which contain information about
several atmospheric constituents
that are each observed in several
spectral elements. The
multiplicity of unknowns and the
redundancy of the data to be handled
lead to the adoption of a retrieval
strategy based upon the
following three choices.
2.4.4.1.3.1.1 Use of Microwindows The redundancy of
information coming from MIPAS
measurements makes it possible to
select a set of narrow (less than 3
cm^{1} width) spectral
intervals containing the best
information on the target
parameters, while the intervals
containing little or no information
can be ignored. The use of selected
spectral intervals, called
'microwindows', allows the
size of analyzed spectral elements
to be limited and avoids the
analysis of spectral regions which
are characterized by uncertain
spectroscopic data, interference by
nontarget species, NonLocal
Thermal Equilibrium (NLTE) and line
mixing effects. More generally,
priority can be given to the
analysis of spectral elements
with most information on the target
species and less affected by
systematic errors, e.g. for VMR retrievals
transitions with weak temperature
dependence can be preferred in order
to minimize mapping of
temperature uncertainties on to
the resulting VMR vertical profiles. 2.4.4.3.
By analyzing the
sensitivity of the radiance to
changes of target parameters, lists
of appropriate microwindows have
been selected for the retrieval
of H_{2}O, O_{3},
HNO_{3}, CH_{4},
N_{2}O and NO_{2} as
well as for the joint retrieval
of pressure and temperature
Clarmann T. v., A. Dudhia,
and CO
Ref. [1.17 ]
. A microwindow
database 2.4.4.4. has been created and is
currently being refined with
respect to minimization of retrieval
errors, following the approach
illustrated in
Clarmann T. v., and G. Echle
Ref. [1.28 ]
.
2.4.4.1.3.1.2 Sequential retrieval of
the species
The unknowns of the retrieval
problem are:
 the observation geometries,
that are identified by the
'tangent
pressures', i.e.
the values of pressure
corresponding to the tangent
points of the limb
measurements (pressure is
the independent
variable and it is used as
the altitude coordinate of
all the profiles) ;
 the profile of temperature;
 the profile of VMR of the
five target species;
 the atmospheric
continuum, that includes
all the emission
sources that are frequency
independent within a microwindow.

zerolevel
calibration correction,
that accounts for
additive microwindow
dependent offsets that could
remain uncorrected after Level
1b 2.4.3.2. processing.
These unknowns are
retrieved using the following
sequence of operations. First
temperature and 'tangent
pressures' are retrieved
simultaneously (p,T retrieval), then
the target species VMR profiles are
retrieved individually in
sequence. The reason of this
approach is that a simultaneous
retrieval of all the species would
require a huge amount of
computer memory, because the size of
the matrices to be handled by the
retrieval is proportional to the
product between the unknown
parameters and the number of
observations. The
feasibility of the simultaneous
retrieval of pressure and
temperature has been investigated in
Carlotti M., and M. Ridolfi
Ref. [1.13 ]
and
Clarmann T. v., A.
Linden,and CO
Ref. [1.19 ]
. Simultaneous p,T retrieval
exploits the hydrostatic equilibrium
assumption, that provides a
relationship between temperature,
pressure and geometrical altitude.
The use of the hydrostatic
equilibrium assumption is discussed
in this section
eq.
2.37 .
The sequence of the target species
retrieval has been determined
according to the degree of their
reciprocal spectral interference and
is: H_{2}O,
O_{3}, HNO_{3},
CH_{4}, N_{2}O and
NO_{2} . In p, T
retrieval the quantities to be
fitted are the 'tangent
pressures', the temperature
profile sampled in correspondence of
the 'tangent pressures'
and the parameters of atmospheric
continuum and zerolevel calibration
correction. In each VMR retrieval, the
fitted quantities are the VMR
altitude distribution of the
considered gas sampled at the
tangent pressures and the parameters
of atmospheric
continuumand instrumental zerolevel
calibration correction.
2.4.4.1.3.1.3 Global Fit analysis of
the LimbScanning sequence
A global fit approach
Carlotti M
Ref. [1.15 ]
is adopted for the
retrieval of each vertical
profile. This means that
spectral data relating to a
complete limb scan sequence are
fitted simultaneously. Compared
to the onionpeeling method
McKee T. B., R. I.
Whitman, J. J. Lambiotte jr
Ref. [1.52 ]
,
Goldman A., R. S. Saunders
Ref. [1.42 ]
, the global fit provides a
more comprehensive exploitation
of the information and a
rigorous determination of the
correlations between atmospheric
parameters at the different
altitudes. Besides, it
permits the full exploitation of
hydrostatic equilibrium
condition and is better
compatible with the modeling of
the finite field of view (FOV) of the instrument.
2.4.4.1.3.2 Mathematics of the retrieval
In the inversion algorithm the
syntetic spectra simulated using a
radiative transfer model (forward
model) through an inhomogeneous
atmosphere are fitted to the
observed spectra. The simulations
are fitted to the observations by
varying the input parameters of the
model (such as pressure,
temperature, VMR, etc...) according
to a nonlinear GaussNewton
procedure. Further details of the
mathematics of the retrieval
model can be found here 2.4.4.1.3.2.1. .
2.4.4.1.3.2.1 Details
The problem of retrieving the
altitude distribution of a
physical or chemical quantity
from limbscanning observations
of the atmosphere falls
within the general class of
problems that requires the
fitting of a theoretical model,
describing the behavior of the
observed system, to a set of
available observations
Ref. [1.65 ]
,
Ref. [1.41 ]
,
Ref. [1.53 ]
,
Ref. [1.61 ]
. The
instrument observes the
radiance emitted by the
atmosphere at different values
of the spectral
frequency and of the
limbviewing angle θ. The
theoretical model (or forward
model ) simulates the
observations through a set of
parameters p
and q:
p represents
the quantities that affect the
radiance but are not
retrieved, and
q the
quantities that are retrieved,
i. e. the distribution profile
of the atmospheric quantity
under investigation.
The retrieval procedure consists
of the search for the set of
values of the parameters
q that produces
the "best"
simulation of the
observations. It has
to be noted that the real
atmospheric profile is a
continuous function, but, since
there will always be a finite
number of measurements
(n), in order to
constrain the problem, the
unknown function is approximated
with a discrete representation
that we indicate with a
vector x of
dimension m, where
m < n. Between
the m discrete points,
an interpolated value is
then used in the forward
model. Assuming a
normal (Gaussian) distribution
for the measurement errors, the
approach generally used for
determining the parameters which
produce the best simulation of
the observation is the
leastsquare fit. This
procedure, deriving from the
theory of the maximum likelihood
estimation
Ref. [1.61 ]
, looks for the
solution x that
minimizes the χ
^{2}
function, defined as the
square summation of the
differences between observations
and simulations, weighted by the
measurement noise.
Given n measurements S_{i}
and the corresponding
simulations F_{i}
(
p
,) calculated by the
forward model using the
assumed profile , and calling
n the
vector of the differences
between observations and
simulations and V
_{n}
the variancecovariance
matrix (VCM) associated to
the vector n,
the quantity to be minimized
with respect to the unknown
parameters x is:
  eq 2.6 
If the
observations are suitably chosen,
matrix V
_{n}
is not singular and its
inverse exists. However, in case
the spectrum is sampled on a grid
finer than 1/(2 * MPD) where MPD is
the Maximum Path Difference, the
inverse of matrix V
_{n}
does not exist. In this case
the generalized inverse
Ref. [1.46 ]
of V
_{n}
is used in equation
eq.
2.6 .
In general, the observations do not
depend linearly on the unknown
parameters x. As a
consequence equation
eq.
2.6 is not a
quadratic function of the unknowns,
and an analytic expression for
the values of the unknowns which
minimizes equation
eq.
2.6 cannot be
determined. However,
sufficiently close to the minimum,
we may assume that the χ
^{2}
function is well approximated
by a quadratic form, obtained
expanding equation
eq.
2.6 in the
Taylor series about the initial
guess profile:
  eq 2.7 
where∇ and
∇
^{2}
indicate respectively the
gradient and the Hessian matrix of
the χ
^{2}
function and is a vector of
dimension m, providing the
correction to be applied to the
assumed value of parameter in order to obtain its
correct value x.
Writing the gradient and
the Hessian matrix of the χ
^{2}
function explicitly, and
indicating with
K the Jacobian
matrix, i.e. a matrix of n
rows and m columns, whose
entry k_{ij}
is the derivative of
simulation i made with
respect to parameter j, we
obtain:
  eq 2.8 
If the problem is
linear or if near the minimum the
residuals have a null average, the
term with / can be neglected in equation
eq.
2.8 and the
value of y
which minimizes the χ
^{2}
function is the GaussNewton
solution
Ref. [1.41 ]
:
  eq 2.9 
Therefore, the
solution matrix of the inverse
problem, defined as the matrix that
calculates the unknowns from the
measured quantities, is equal
to:
  eq 2.10 
If the hypothesis
of linearity is not satisfied, with
this procedure the minimum of the
χ
^{2}
function is not reached
but only a step is done toward the
minimum. The vector
, with
y computed using equation
eq.
2.9 ,
represents only a better estimate of
the parameters than . In this case the whole
procedure must be reiterated
starting from the new estimate of
the parameters (Newtonian iteration)
and equation
eq.
2.9 has to be
written as:
  eq 2.11 
where
iter indicates the
iteration index, the result of the
previous iteration,
the Jacobian relative
to the profile , the residuals.
Convergence criteria are therefore
needed in order to establish when a
value close enough to the minimum of
the χ
^{2}
function has been reached.
However, this procedure is
successful only in the case of
sufficiently weak
nonlinearities. If we start from a
position of the χ
^{2}
function that is far from the
minimum, its second order
expansion may be a poor
approximation of the shape of the
function, so that the calculated
correction can be misleading, and
increase rather than decrease
the residuals. For this reason, a
modification of the GaussNewton
method, the Levenberg
Ref. [1.47 ]
 Marquardt
Ref. [1.51 ]
,
Ref. [1.56 ]
method, is used.
The modification involves the
introduction in equation
eq.
2.11 of a
factor λ which reduces the
amplitude of the parameter
correction vector.
  eq 2.12 
The factor λ
is initialized to a userdefined
(less than 1) number, and during the
retrieval iterations it is increased
or decreased depending on
whether the χ
^{2}
function increases or
decreases. At convergence the
Levenberg  Marquardt method
provides the same solution that is
obtained with equation
eq.
2.11
because they aim at the same minimum
of the χ
^{2}
function. If V
_{n}
is a correct estimate of the
errors of the observations and the
real minimum of the χ
^{2}
function is found, the
quantity defined by equation
eq.
2.6 has an
expectation value of (n 
m) and a standard deviation
equal to . The value of the
quantity provides therefore a
good estimate of the quality of the
retrieval. Values of that deviate too much
from unity indicate the presence of
incorrect assumptions in the
retrieval. The errors
associated with the solution of the
inversion procedure can be
characterized by the
variancecovariance matrix V
_{x}
of x given
by:
  eq 2.13 
where D
_{c}
and K
_{c}
are respectively the solution
matrix and the Jacobian matrix
evaluated at convergence.
Matrix V_{x}
maps the experimental
random errors onto the uncertainty
of the values of the retrieved
parameters. Actually, the square
root of the diagonal elements of V_{x}
measures the root mean
square (r.m.s.) error of the
corresponding parameter. The
offdiagonal element (V
_{x}
)_{ij}
of matrix V_{x}
, normalized to the square
root of the product of the two
diagonal elements (V
_{x}
)_{ii}
and (V
_{x}
)_{jj}
, provides the correlation
coefficient between parameters
i and j.
According to the theory of maximum
likelihood, if an apriori
information on the unknown
parameters is available, (x
_{a} and V
^{1}
_{a} being respectively the
vector containing the apriori
values of the unknown and its VCM) the χ
^{2}
expression to be minimized is
Ref. [1.59 ]
:
  eq 2.14 
and equation
eq.
2.11 becomes:
  eq 2.15 
while equation
eq.
2.13 becomes:
  eq 2.16 
If we assume that
the complementary information
consists of a vector n_{1}
with a VCM
connected by the Jacobian K_{1}
to the unknowns
y, equation
eq.
2.15 can be
written as:
  eq 2.17 
This expression,
which is implicitly used in equation
eq.
2.9 for
the combination of the information
provided by independent microwindows, is
used for the exploitation of
nonradiometric information
(external information). In section 5.1. the
problems associated with the use of
external information on the unknown
parameters will be reviewed on
the light of the choices implemented
in the MIPAS Level 2 retrieval
modules.
2.4.4.1.3.3 Main components of the retrieval
As already stated in the previous
section, the objective of the
retrieval program is the
determination of the atmospheric
parameters that better fit the
simulations to the
observations. Starting
from some firstguess values
of the unknown parameters 2.4.4.5.
and using information on
observation geometry and
instrumental characteristics, the
forward model computes the simulated
spectra, which are compared with
the measured spectra provided by
MIPAS Level 1b
processor 2.4.3.2. . The difference
between simulated and measured
spectra provides the vector of the
residuals used to evaluate the cost
function ( χ;
^{2}
) to be minimized (see mathematical
details here) 2.4.4.1.3.2.1. . A new
profile is generated by modifying
the first guess with the corrections
provided by the GaussNewton
formula. The inversion requires
the knowledge of the Jacobian
matrix. The improved profile can be
used as new guess for generating
simulated spectra which are again
compared with the measured ones.
The iterative procedure stops when
the convergence criteria are
fulfilled. The main
components of the retrieval
algorithm are therefore:
2.4.4.1.3.3.1 Forward model
The purpose of the forward model
is to simulate the spectra
measured by the instrument in
the case of known atmospheric
composition. The signal
measured by the spectrometer is
equal to the atmospheric
radiance which reaches the
spectrometer modified by the
instrumental effects, mainly
due to the finite spectral
resolution and the finite Field
of View (FOV) of the
instrument. The atmospheric
radiance that reaches the
instrument when pointing to the
limb at tangent altitude z_{t}
is calculated by means of
the radiative transfer equation:
  eq 2.18 
where x
is the position along the line of
sight between the observation point x_{0}
and the point x_{i}
at the farthest extent of the
line of sight, B(, x) is the
source function, c(,x) is the
absorption crosssection,
η(x) is the number
density of absorbing molecules.
The exponential term represents the
atmospheric transmittance between
x and x_{0}
. In the case of local
thermodynamic equilibrium B(, x) is
equal to the Planck function.
The measured signal S(,z_{t})
is simulated convolving the
atmospheric limb radiance with the Apodized
Instrument Line Shape
(AILS, defined in section 2.4.4.1.3.3.1.5. ) and with
the MIPAS FOV function (FOV(z_{t}
), discussed in section 2.4.4.1.3.3.1.6. ):
  eq 2.19 
The computation of
equation
eq.
2.18 and equation
eq.
2.19 requires
many operations that must be
repeated for several variables,
each with numerous values. The
search for a sequence of operations
that avoids repetition of the same
calculations and minimizes the
number of memorized quantities is
the first objective of the
optimization process. The
following sequence of operations has
been chosen:

 1. definition of the
frequency grid in
which the
atmospheric radiance
is calculated (section 2.4.4.1.3.3.1.1. );

2. definition of
the tangent
altitude grid of
the simulations
and of an
atmospheric
layering (common
to all the
simulations)
used
for discretising
their radiative
transfer
integral;
determination of
the
'paths'
(i.e.
combination of
layer and limb
view) that
require a
customized
calculation of
the values of
pressure and
temperature and
definition of
the values of
pressure
and temperature
characteristics
of each of these
'paths'
( section 2.4.4.1.3.3.1.2. );
3. computation
of absorption
crosssections,
relating to all
the selected
paths and all
the gases ( section 2.4.4.1.3.3.1.3. );
The implementation of each
of these operations is described
below.
2.4.4.1.3.3.1.1 Dependence on
spectral frequency
Limb radiance spectra contain
spectral features varying
from the sharp, isolated,
Dopplerbroadened lines at
high altitudes to wide,
overlapping,
Lorentzbroadened lines at
low altitudes. In order to
resolve the sharp features
at high altitude, a
'fine grid' of
spectral resolution of the
order of 0.0005
cm^{1} is required.
Furthermore the overlapping
wings at low altitudes
require the grid to be
ubiquitous. The choice of a
small spacing implies a
large number of spectral
points and an equally large
number of
calculations.
However, even if a fine grid
of 0.0005 cm^{1} is
needed, not all the points
of this grid are equally
important for the
reconstruction of the
spectral distribution and
the full radiative transfer
calculation needs only be
performed for a subset of
fine grid points, the
remaining spectrum being
recovered by interpolation.
This subset is denoted the
'irregular grid'
and depends on the microwindow
boundaries and on the
Instrument Line Shape.
In principle it is possible
to determine an optimal
irregular grid for each
tangent altitude, taking
advantage of the particular
line broadening
associated with each
observation geometry.
However, while this might
lead to a reduction in total
number of fine grid points
for which radiative
transfer calculations are
required, this prevents the
savings obtained with
calculations common to all
tangent altitudes. As a
consequence, an irregular
grid is determined which is
valid for all tangent
altitudes, and therefore
common to all observation
geometries. The
number of points that are
selected for the irregular
grid depends upon the
interpolation law that is
used for the reconstruction
of the spectral
distribution. Several
interpolation functions were
investigated, but since the
same interpolation is also
required for the
calculation of Jacobian,
which, unlike the radiances,
can have negative values,
logarithmic and inverse
interpolations had to be
discarded and finally a
simple linear interpolation
was chosen. The
procedure for the generation
of the 'irregular
grid' is to start with
a set of complete fine grid
limb radiance spectra for
different tangent altitudes.
Then, for each point, a
`cost' function is
determined, representing the
maximum radiance error (at
any altitude, after AILS
convolution) that would
arise if the point were to
be replaced by an
interpolated value. The
point with the smallest
interpolation cost is then
eliminated and the process
repeated for the remaining
grid. The iteration
continues until no further
points can be removed
without exceeding a maximum
error criterion, chosen to
be 10 % of the
noiseequivalent signal
radiance. In figure2.23
an example of the high
resolution radiance of a
representative microwindow
computed for the irregular
grid points, as well as its
convolved radiance and the
difference between the
original and interpolated
radiance, is shown.
Typically it is found that
only 510 % of the complete
fine grid is sufficient for
a satisfactory
reconstruction of the
spectral distribution.

Figure 2.23 44 km tangent height spectral radiance from a microwindow (frequency range 693.45  693.725 cm<sup>1</sup>) selected for p, T retrieval. The upper plot shows the highresolution radiance (solid line), the irregular grid points (+), and the convolved radiance (dashed line). The lower plot shows the difference between the original and interpolated radiances, on both the highresolution grid (solid line, left axis) and the convolved radiances (dashed line, right axis). 
2.4.4.1.3.3.1.2 Raytracing and
definition of the
atmospheric layering
The radiative transfer
integral, given by equation,
is a path integral along the
line of sight in the
atmosphere. The line of
sight is determined by the
viewing direction of the
instrument and, due to
refraction, is not a
straight line, but bends
towards the earth. The
refractive index of air is a
function of both pressure
and temperature (the
dependence on frequency is
negligible in the spectral
region of MIPAS
measurements, as well as
the dependence on water
vapor, since MIPAS does not
penetrate in the lower
troposphere) and can be
determined with the Edlen
Edlen B.
Ref. [1.32 ]
model. Since the Earth
is assumed locally
spherical, the local radius
of curvature being
determined by the simplified
WGS84 model
Department of Defense
Ref. [1.25 ]
which has been
adopted for the Earth shape,
the atmospheric layering is
spherical too. In these
conditions, the optical path
x is linked to
the altitude r by
the following expression:
  eq 2.20 
with
r altitude referred to
the earth center, n(r)
refractive index and r_{t}
tangent altitude.
Equation
eq.
2.20 has a
singularity at the tangent point
(r = r_{t}
), however the singularity
can be removed by changing the
integration variable from
r to. It
results that:
  eq 2.21 
and the limit
of this expression for
r > r_{t}
can be computed
analytically considering the
dependence of the refractive
index, the pressure and the
temperature on the altitude.
The path integral equation is
computed as a summation over
a set of discrete layers. A
common set of layers is defined
for all the spectra of the
sequence. The boundaries of
these layers are defined in
correspondence of the grid of
the simulated tangent altitudes.
We recall that spectra are
simulated at the tangent
altitudes of the limb scan
sequence and at some additional
tangent altitudes that are used
for the FOV
convolution (as discussed in section 2.4.4.1.3.3.1.6. ).
Extra levels are added to the
simulation grid as long as
either the variation of
temperature or the Voigt
halfwidth of a reference line
across each layer are greater
than a maximum threshold
supplied by the user.
For each layer that results from
this process, appropriate
'equivalent' pressure
and temperature, namely the
CurtisGodson
Houghton J. T.
Ref. [1.45 ]
quantities, have to be
determined. These are calculated
weighting the pressure and
temperature along the raypath
with the number density of each
absorbing gas. This technique
allows a coarse discretisation
of the atmosphere to be
implemented. The equivalent
value of parameter G
(pressure or temperature
respectively) relative to the
gas g, the layer
l and the limb view
t, is calculated as:
  eq 2.22 
where
z is the
altitude, and are the heights of
the boundaries of the
layer, is the VMR of the
gth gas, x
_{t} is the line of
sight characterized by the
tangent altitude z_{t}
, is the air number
density, and is the slant
column relative to the
considered gas, layer and limb
view:
  eq 2.23 
It has to be
noted that in principle,
CurtisGodson pressures and
temperatures have to be computed
for each gas, each layer and
each limb view. However,
when the approximations of flat
layers and straight line of
sight are valid, both the
numerator and the denominator of
equation
eq.
2.22 are
proportional to the secant of
the angle θ
between the line of sight and
the vertical direction and the
same values of p_{e}
and T_{e}
are obtained independently
of the limb angle. We have
verified that the use of secant
law approximation causes very
small errors at all altitudes,
except at the tangent layer, as
it is shown in figure2.24 .
In this approximation, since a
layering common to all spectra
of the scan is defined, it is
sufficient to calculate the
values of p_{e}
and T_{e}
for all the layers of the
lowest limb view, and only for
the lowest layer of the other
limbviews, all the other layers
being characterized by the same
values of p_{e}
and T_{e}
as the lowest limb
view. This
optimization is crucial, not
only because less equivalent
pressures and temperatures have
to be calculated, but mainly
because absorption
crosssections have to be
calculated and stored for a
smaller number of p, T
combinations.

Figure 2.24 Differences between CurtisGodson pressures and temperatures of the layers of the limb view with tangent altitude 30 km and the corresponding quantities of the same layers in the case of the lowest limb view (tangent altitude 6 km). On the left axis the percent differences on equivalent pressures (squares) are reported, on the right axis the absolute differences on equivalent temperature (triangle) are shown. The layers are counted from the tangent layer at 30 km and are 2 km thick. 
