## Abstract

The trend of stratospheric water vapor as a function of latitude is estimated by the MIPAS measurements by means of a new method that uses the measurement space solution. The method uses all the information provided by the observations avoiding the artifacts introduced by the a priori information and by the interpolation to different vertical grids. The analysis provides very precise values of the trends that, however, are limited by a relatively large systematic error induced by the radiometric calibration error of the instrument. The results show in the five years from 2005 to 2009 a dependence on latitude of the stratospheric (from 37 to 53 km) water vapor trend with a positive value of (0.41 ± 0.16)%yr^{−1} in the northern hemisphere and less than 0.16%yr^{−1} in the southern hemisphere.

© 2011 OSA

## 1. Introduction

Water vapor is the most important gaseous source of infrared opacity in the atmosphere, it accounts for about 60% of the natural greenhouse effect for the clear skies [1], and provides the largest positive feedback in model projections of climate change [2]. Therefore, water vapor variability is an important issue in the discussion of global climate change [3] and in particular the variability of stratospheric water vapor has important radiative and chemical consequences that impact the global surface climate change [4].

An increase of roughly 1% per year in stratospheric water vapor content has been observed during the last half of the 20th century [5, 6], with a more convincingly documented increase during the 1980s and most of the 1990s than earlier. However, an updated trend analysis [7] of water vapor in the lower mid-latitude stratosphere from Boulder balloon measurements and from HALOE (Halogen Occultation Experiment) [8] spaceborne observations provides trend estimates for the period 1980-2000 that are up to 40% lower than previously reported.

Methane oxidation is a major source of water in stratosphere, and has been increasing over the industrial period, however, the observed trend in stratospheric water vapor during the last half of the 20th century is too large to be attributed to methane oxidation alone [5, 9].

The temperatures near the tropical tropopause should control the stratospheric water vapor content according to the equilibrium thermodynamics, importing more water vapor into the stratosphere when temperatures are warmer. However, tropical tropopause temperatures have cooled slightly over the period of the stratospheric water vapor increase [10, 11]. Other mechanisms have been proposed to explain the increase of the stratospheric water vapor occurred in the second half of 20th century, but so far the driving causes of this increase are unknown.

The upward trend of stratospheric water vapor decreased in the last half of the 1990s with a near-zero trend between 1996 and 2000 [12, 13]. Furthermore, at the end of 2000 there was a dramatic drop of about 10% of stratospheric water vapor [13]. The trend analysis reported in [14] extends until spring 2008 and it shows that a minimum was approximately reached between 2004 and 2006 and an increase is observed afterwards.

The drop in stratospheric water vapor that occurred at the end of 2000 is thought to have slowed the rate of increase in global surface temperature over 2000-2009 by about 25% compared to that which would have occurred due only to carbon dioxide and other greenhouse gases [4]. On the other hand the increase in stratospheric water vapor occurred between 1980 and 2000 would have enhanced the decadal rate of surface warming during the 1990s by about 30%. These considerations show that stratospheric water vapor is an important driver of decadal global surface climate change and the accurate determination of its variability is, therefore, essential for interpreting global changes and for making reliable future projections.

Trend determination must avoid biases due to the systematic errors of the measurement technique and to time and geographical natural variability.

In this paper we make an estimation of the trend of stratospheric water vapor as a function of latitude in the five years from 2005 to 2009 using an homogeneous set of measurements acquired by the Michelson Interferometer for Passive Atmospheric Sounding (MIPAS) onboard the ENVISAT satellite [15], and discuss all the precautions taken to either avoid or correct for possible biases. Five years is a short time for the determination of a trend, but it provides a valuable example of the data processing that could be used for an accurate determination of the variability of the atmospheric parameters in view of an early detection of the on going processes.

The use of measurements performed by a satellite instrument to derive the trend values has the advantage to exploit a set of internally consistent measurements with global coverage for several years. However, since the remote sensing measurements do not directly observe the quantity of which we estimate the trend, the procedure used for the retrieval must be critically analyzed and some attention has to be paid to the analysis of the quantities retrieved from the observations. Indeed, when the vertical profile of an atmospheric constituent is retrieved from remote sensing measurements the retrieved profile is related to the true profile by a relationship that depends on the a priori information (including regularization constraints) [16] used to make the retrieval well-conditioned. The relationship between the retrieved and the true profile is described by means of the averaging kernels, which, if they change from measurement to measurement, can introduce an artificial variability in the data and affect the trend determination. Indeed, if the averaging kernels change it means that we are analyzing the time dependence of a parameter that does not represent always the same physical quantity, and the variability that we observe can be due to the fact that we are measuring different quantities at different times.

The use of retrieved profiles to estimate trends poses also the problem that the vertical retrieval grid can be different for different measurements and, therefore, an interpolation to a common grid of all the measurements is needed. This interpolation to a common vertical grid determines always a loss of information [17] that degrades the quality of the data.

In order to avoid the problems posed by the a priori information and by the interpolation on a common grid, we use the approach of the *measurement space solution* (MSS) [18–20], which also makes sure that all the information provided by the observations is exploited. Indeed the MSS approach has the advantage that the solution can be represented on a vertical grid as fine as desirable without relying on any a priori information.

Using the MSS approach we estimate the stratospheric water vapor trend and the seasonal variability as a function of latitude and for all the globe with an analysis of all the error components.

In Section 2 we describe the mathematical formalism for the estimation of the trends. Firstly we discuss the method that should be adopted when using the retrieved vertical profiles of water vapor and then we describe the MSS approach that is used in this paper. In Section 3 we apply the MSS approach to the MIPAS measurements and describe the results. Finally in Section 4 we draw the conclusions.

## 2. Theory

#### 2.1. Definition of the problem

We represent the vertical profile of an atmospheric variable with a column vector **x** of *n* elements corresponding to a vertical (in altitude or pressure) grid fine enough to represent adequately the altitude dependence of the variable. In this paper we consider the estimation of trends of concentrations of water vapor and, therefore, the elements of the profile **x** represent volume mixing ratios (VMRs) of water vapor, however, the same method can be used to estimate trends of any atmospheric variable. In general, the variable *c*(*t*) (depending on time) of which we wish to estimate the trend is a linear function of the vertical profile **x** which can be represented by means of a dot product between a vector **w** and **x**:

*denotes the transposition operation. The vector*

^{T}**w**defines the vertical range to which the parameter

*c*(

*t*) refers and generally is a peaked function of altitude satisfying the normalization condition In this case we can interpret

*c*(

*t*) as the value of VMR at the altitude of the maximum of

**w**with a vertical resolution given by the width of the peak. Larger widths of the peak correspond to smaller measurement errors of the parameter

*c*(

*t*), therefore, the width of the peak has to be selected on the basis of the trade off between vertical resolution and precision.

In case that we are interested in the vertical column, we take the *i ^{th}* element of

**w**,

*w*, equal to the air column in the

_{i}*i*layer.

^{th}The problem that we face is the estimation of the parameter *c*(*t*) as a function of time.

#### 2.2. Estimation of the parameter c(t) from a retrieved profile

In this subsection we establish the formalism for the estimation of the parameter c(t) from a retrieved profile. While this approach is not used in this paper, it is useful to explain, from a theoretical point of view, why the MSS approach provides better performances. The reader not interested in the formalism can skip this subsection.

The parameter *c*(*t*) defined in Eq. (1) can be exactly calculated when we know the true profile **x**(*t*). In the real case we have a profile **x̂**(*t*) retrieved from a measurement which is related to **x**(*t*) by the equation [16]:

*ε*(

*t*) contains the errors on the retrieved profile propagated from the errors that affect the observations,

**A**is the averaging kernel matrix,

**x**

_{0}(

*t*) is a profile close enough to the true profile

**x**(

*t*) in such a way that the linear approximation of the averaging kernel is appropriate and

**x̂**

_{0}(

*t*) is the retrieved profile when

**x**

_{0}(

*t*) is the true profile and the observations are not affected by errors.

**x̂**

_{0}(

*t*) can be obtained with a simulated retrieval from error-free observations computed with the forward model assuming the

**x**

_{0}(

*t*) profile as the true profile. The vertical grid on which

**x̂**(

*t*) and

**x̂**

_{0}(

*t*) are represented is the retrieval grid which may depend on the measurement and retrieval technique and, therefore, is in general different from the vertical grid on which

**x**(

*t*) and

**x**

_{0}(

*t*) are represented.

The quantity *c*(*t*) can be estimated substituting in Eq. (1) the retrieved profile **x̂**(*t*) in place of the true profile **x**(*t*) after that an interpolation is made from the retrieval grid to the grid of **x**(*t*) with a suitable interpolation matrix **H**. Proceeding in this way and using Eq. (3) we obtain:

*ĉ*(

*t*) instead of

*c*(

*t*) is

**w**

*[*

^{T}**Hx̂**

_{0}(

*t*) –

**x**

_{0}(

*t*)] can be calculated and used to correct

*ĉ*(

*t*). Consequently our estimation of

*c*(

*t*) is affected by two errors:

**w**

*(*

^{T}**HA**–

**I**)[

**x**(

*t*) –

**x**

_{0}(

*t*)], that is due to the smoothing error [16] of the measurement method, and

**w**

^{T}**H**

*ε*(

*t*), that is due to the errors of the observations.

The smoothing error component can be estimated as suggested in [16]. The statistics of this error can be calculated from the mean **x**
* _{e}* and the covariance matrix

**S**

*of an ensemble of states that comprehensively describes both the time and the geographical variability of the profile. The mean of this error is*

_{e}**w**

*(*

^{T}**HA**–

**I**)[

**x**

*(*

_{e}*t*) –

**x**

_{0}(

*t*)] (which gives a bias term) and the variance about

**x**

*is*

_{e}**w**

*T*(

**HA – I**)

**S**

*(*

_{e}**HA – I**)

^{T}**w**. As stated in [16] the correct estimate of the smoothing error contribution requires to know the actual statistics (

**x**

*and*

_{e}**S**

*) of the fine structure of the profile. If this is not known and we use some*

_{e}*ad hoc*matrix that has been constructed as a reasonable a priori constraint in the retrieval we can make a significant error in the estimation of this contribution. The smoothing error is an elusive and often underestimated contribution which, as it will be later discussed, is due to those features of the observed quantity that are not measured by the data analysis process.

We can try to reduce the contribution due to the smoothing error replacing **w**
^{T}**H** in Eq. (4), which is the transformation used to determine our parameter *ĉ*(*t*) from the retrieved profile **x̂**(*t*), with an optimized transformation **v**
* ^{T}*, where

**v**is a vector (with a number of elements equal to the number of points of the retrieval grid) to be determined:

In this case, the error that we make estimating the quantity *ĉ*(*t*) instead of *c*(*t*) is

**v**

^{T}**A – w**

*)[*

^{T}**x**(

*t*) –

**x**

_{0}(

*t*)] and it can be minimized choosing the vector

**v**that minimizes the norm of the vector (

**A**

^{T}**v – w**). The vector

**A**

^{T}**v**is the linear combination of the rows of the matrix

**A**(the averaging kernels) with the coefficients given by the elements of

**v**. This property suggests to write the vector

**w**as the sum of the component

**w**

*of*

_{A}**w**in the space spanned by the rows of

**A**(that we indicate with {

*A*}) plus the component

**w**

_{A}_{⊥}of

**w**in the orthogonal complement space in ℝ

*(including all the vectors made of ordered*

^{n}*n*-tuples of real numbers) of {

*A*} (that we indicate with {

*A*

^{⊥}}):

**A**

^{T}**v – w**

*) and*

_{A}**w**

_{A}_{⊥}are orthogonal, the square norm of (

**A**

^{T}**v – w**) is given by

**A**

^{T}**v – w**)| is obtained choosing the vector

**v**=

**v**

*for which*

_{min}**A**

^{T}**v**

*=*

_{min}**w**

*and Eq. (9) becomes:*

_{A}**v**=

**v**

*the error contribution due to the smoothing error is proportional to the component of*

_{min}**w**that belongs to the space {

*A*

^{⊥}} and, therefore, is zero in the case that the vector

**w**belongs to the space {

*A*} (that is if

**w**is a linear combination of the averaging kernels). If

**w**does not belong to {

*A*} the statistics of this error can be estimated (as described above) with the mean value ${\mathbf{\text{w}}}_{A\perp}^{T}[{\mathbf{\text{x}}}_{0}(t)\hspace{0.17em}-\hspace{0.17em}{\mathbf{\text{x}}}_{e}(t)]$ (which gives a bias term) and the variance ${\mathbf{\text{w}}}_{A\perp}^{T}{\mathbf{\text{S}}}_{e}{\mathbf{\text{w}}}_{A\perp}$. This means that the smoothing error is an “incompleteness error” due to the fact that the space generated by the averaging kernels {

*A*} does not coincide with the complete space ℝ

*. Indeed when {*

^{n}*A*} coincides with ℝ

^{n}**w**

_{A}_{⊥}is the null vector and the smoothing error is zero.

Once the vector **v**
* _{min}* has been determined, the contribution
${\mathbf{\text{v}}}_{\mathit{\text{min}}}^{T}{\widehat{\mathbf{\text{x}}}}_{0}(t)\hspace{0.17em}-\hspace{0.17em}{\mathbf{\text{w}}}^{T}{\mathbf{\text{x}}}_{0}(t)$ can be calculated and used to correct

*ĉ*(

*t*) and the contribution ${\mathbf{\text{v}}}_{\mathit{\text{min}}}^{T}\varepsilon (t)$, due to the errors on the observations, can be represented by the variance ${\mathbf{\text{v}}}_{\mathit{\text{min}}}^{T}\mathbf{\text{S}}{\mathbf{\text{v}}}_{\mathit{\text{min}}}$, where

**S**=<

*ε*(

*t*)

*ε*(

*t*)

*> is the covariance matrix of*

^{T}**x̂**(

*t*) (the symbol < ... > representing the mean value).

The averaging kernel matrix is given by [16]:

where**G**is the gain matrix (which describes the sensitivity of the retrieval to the observations) and

**K**is the Jacobian matrix of the forward model calculated in

**x**

_{0}(

*t*).

Equation (11) shows that the rows of **A** (the averaging kernels) are linear combinations of the rows of **K** with the coefficients given by the elements of the rows of **G**. Therefore, the space {*A*} is a subspace of the space spanned by the rows of **K** [21] which is referred to as the *measurement space* [18–20] and that we indicate with {*K*}. Consequently, the incompleteness error (smoothing error) is smaller if we use {*K*} in place of {*A*}. This statement can be deduced from the following consideration. Equation (3) shows that the true profile contributes to the retrieved profile by means of its projections on the rows of **A**, that is by way of its component in the space {*A*}. Therefore, the retrieved profile **x̂**(*t*) allows to retrieve the component of **x**(*t*) in the space {*A*}. However, the observations provide the projections of the true profile **x**(*t*) on the rows of **K** [18] and, therefore, they contain the information for retrieving the component of **x**(*t*) in the space {*K*} which includes the space {*A*} and is more complete. Consequently the use of the retrieved profile **x̂**(*t*) (that provides the solution in {*A*}) to estimate *c*(*t*) may not exploit all the information contained in the observations (that provide the solution in {*K*}).

Another problem that arises with the use of the retrieved profile to estimate trends is that generally the averaging kernels are given to the user as square matrices, i.e. they refer to the retrieval grid both for the rows and for the columns. If the retrieval grid is not adequate for the representation of the real profile, some hypothesis has to be made to transform the averaging kernels from the retrieval grid into the grid on which the true profile is represented. The approximations made in this hypothesis may introduce a degradation in the information content of the data.

The optimal way to estimate *c*(*t*) exploiting all the information contained in the observations (with a minimum incompleteness error) and without any problem due to interpolation between different vertical grids is to use the MSS method [18]. Indeed the MSS provides the solution in the space {*K*} and on a vertical grid as fine as desirable.

#### 2.3. Estimation of the parameter c(t) from the MSS

In the MSS method the vertical grid of the retrieved profile and of the Jacobian matrix **K** is not constrained, therefore, we can use a fine vertical grid of *n* points equal to that used to represent the profile **x**(*t*). Since the ℝ* ^{n}* space can be split into the direct sum of the measurement space {

*K*} and of its orthogonal complement {

*K*

^{⊥}} (referred to as the

*null space*[18]), we can write:

**x**

*(*

_{K}*t*) and

**x**

_{K}_{⊥}(

*t*) belong to the spaces {

*K*} and {

*K*

^{⊥}}, respectively.

**x**

*(*

_{K}*t*) and

**x**

_{K}_{⊥}(

*t*) can be expressed as: where

**V**is a matrix whose columns are an orthonormal basis of {

*K*},

**W**is a matrix whose columns are an orthonormal basis of {

*K*

^{⊥}} and

**a**(

*t*) and

**b**(

*t*) are the projections of

**x**(

*t*) on these orthornormal bases: The component

**x**

*(*

_{K}*t*), expressed in terms of

**V**and

**a**(

*t*), is the MSS and can be derived from the observations following the procedure described in [18, 20]. The elements of the vector

**a**(

*t*) are characterized by noise errors (propagated from the noise on the observations) which are uncorrelated with each other (the covariance matrix of

**a**(

*t*) is a diagonal matrix by construction) and these errors are used to sort the elements of

**a**(

*t*): the error increases with the index of the element. The elements of

**a**(

*t*) can be affected also by some systematic errors due to the uncertainty on some instrumental parameters. The analysis of these systematic errors and of how they affect the estimation of the trends is discussed in subsection 3.4.

Substituting Eq. (12) into Eq. (1) and writing the vector **w** as the sum of the component **w**
* _{K}* in the space {

*K*} plus the component

**w**

_{K}_{⊥}in the space in {

*K*

^{⊥}} we obtain:

**x**

_{K}_{⊥}(

*t*), we can only obtain an estimation of the term ${\mathbf{\text{w}}}_{K\perp}^{T}{\mathbf{\text{x}}}_{K\perp}(t)$ using some a priori information. If we choose the vector

**w**belonging to the measurement space then

**w**

*=*

_{K}**w**,

**w**

_{K}_{⊥}= 0 and using the MSS we can estimate exactly

*c*(

*t*). Since the measurement space depends on the Jacobian matrix which is in general different for different observations, it is not always possible to choose a vector

**w**that belongs to all the measurement spaces of the different observations. We choose a vector

**w**which almost completely lies in all the measurements spaces, that corresponds to the condition |

**w**

_{K}_{⊥}| << |

**w**

*| ∼ |*

_{K}**w**|, and, using a climatological profile, estimate the error that we make neglecting the term ${\mathbf{\text{w}}}_{K\perp}^{T}{\mathbf{\text{x}}}_{K\perp}(t)$. In agreement with what said in subsection 2.2, we refer to this error as the ”incompleteness error” because it arises from the fact that the measurement space does not coincide with the complete space ℝ

*.*

^{n}The procedure has to adopt a strategy to select the number of components to represent the MSS. Indeed the components of the MSS (given by the columns of **V** weighted with the elements of **a**(*t*)) are affected by noise errors, therefore, increasing the number of components of the MSS we increase the noise error on the term
${\mathbf{\text{w}}}_{K}^{T}{\mathbf{\text{x}}}_{K}(t)$. On the other hand, if we represent the MSS with a small number of components the contribution of the term
${\mathbf{\text{w}}}_{K\perp}^{T}{\mathbf{\text{x}}}_{K\perp}(t)$ and the incompleteness error increase. Therefore, the selection of the number of components to represent the MSS is based on the identification of the best compromise between these two competing requirements.

## 3. Estimation of water vapor trends from MIPAS measurements

#### 3.1. MIPAS measurements

We applied the procedure described in subsection 2.3 to a set of measurements acquired by the MIPAS instrument in the period 2005-2009. MIPAS [15] is a Fourier-transform spectrometer operating in the middle infrared that observes the atmospheric emission at the limb for the retrieval of the vertical profiles of several minor atmospheric constituents. The MIPAS measurements are acquired in the nominal measurement mode adopted after January 2005, with an unapodized spectral resolution of 0.0625 cm^{−1} (corresponding to a maximum optical path difference equal to 8 cm) at 27 tangent heights covering the altitude range from 7 to 72 km. The code used for the determination of the linearization point (needed for the MSS calculation [18, 20]) is that adopted by the European Space Agency (ESA) for the operational retrieval [22–25]. It uses a non-linear least-squares fit of the observed spectra with forward model simulations to retrieve the vertical profiles of pressure, temperature, water vapor, ozone, nitric acid, methane, nitrous oxide and nitrogen dioxide. The microwindow approach, described in [26], is adopted and of the 27 spectra only a subset of spectral points containing the relevant information for water vapor retrieval is used.

The results described in this paper refer to 265448 measurements acquired in the five years from January 2005 to December 2009 and made available by ESA for validation purposes. A larger number of measurements is available in the three years from 2005 to 2007 (a mean value of 6350 measurements per month) with respect to the years 2008 and 2009 (a mean value of 1530 measurements per month). No measurements were acquired in the nominal mode in September, October and November 2005 and in February and April 2006, therefore, these months are not included in our analysis.

#### 3.2. Estimation of the parameter c(t)

The vertical grid on which we calculated the MSS is the ECMWF (European Centre for Medium-Range Weather Forecasts) pressure grid with 91 levels for a surface pressure of 1013.250 hPa. It extends from 1013.250 to 0.02 hPa and its detailed description can be found in [27].

For each measurement we performed the retrieval with the code developed for the operational retrieval and interpolated the retrieved profile to the 91 levels pressure grid. This interpolated profile was used as linearization point to calculate the MSS.

The vector **w** was chosen with the objective of selecting a component of the water vapor profile that is extensively determined (with minimum incompleteness error) in each measurement space of all the analyzed measurements and, at the same time, is little affected by the natural variability of the profile distribution. Since the observation geometry of MIPAS allows accurate measurements in stratosphere [23], we chose a Gaussian function of the logarithm of the pressure peaking at 2 hPa and extending over all the stratosphere:

*N*is a normalization term that we determine using Eq. (2),

*p̄*= 2 hPa and

*σ*= 1.6*log(2). In the altitude domain

**w**is approximately a gaussian centered at 45 km with a standard deviation of 8 km. Figure 1 shows the vector

**w**as a function of pressure. We verified that for 97% of the available measurements |

**w**

*| (calculated with the number of components that minimizes the quadratic sum of the noise and the incompleteness errors) is larger than 99% of |*

_{K}**w**|, showing that the component of the water vapor profile along

**w**is almost completely measured in most of the considered measurements.

For each measurement we calculated the term ${\mathbf{\text{w}}}_{K}^{T}{\mathbf{\text{x}}}_{K}(t)$ (see Eq. (17)) with its noise error and estimated the term ${\mathbf{\text{w}}}_{K\perp}^{T}{\mathbf{\text{x}}}_{K\perp}(t)$ (incompleteness error) using a climatological profile [28] corresponding to the season and the latitude of the measurement. For each measurement we selected the number of components which minimizes the quadratic sum of the noise and the incompleteness errors. In the analyzed measurements the number of selected components has a mean value of 20.5 and a standard deviation of 4.4. These values compare with a mean value and a standard deviation of the number of measurements per limb scan equal to 21.1 and 5.4, respectively.

In Fig. 2 we report the histogram of the percentage incompleteness errors obtained from the 265448 analyzed measurements with a bin width of 0.05%. This histogram shows that the incompleteness error is made of a random component superimposed on a bias term that we estimated by the mean value. We also calculated the histogram reported in Fig. 2 separately for each year from 2005 to 2009 verifying that no change with time occurred in the incompleteness error statistics in the analyzed time period.

In order to consider only data measured with large precision we have filtered out the data for which either the noise error is larger than 2% or the incompleteness error is larger than 0.6% (this threshold allowed to reduce the bias term to less than 0.1%) or |**w**
_{K}_{⊥}| is larger than two parts per thousand of |**w**|. The remaining data consist in 184396 measurements (69.47% of the all analyzed measurements) for which *c*(*t*) was estimated as equal to
${\mathbf{\text{w}}}_{K}^{T}{\mathbf{\text{x}}}_{K}(t)$.

#### 3.3. Estimation of the trends

We performed an analysis of the trend for latitude bands. We divided the globe in 17 latitude bands: 15 latitude bands from –75° to +75° with 10° steps and two polar bands from –90° to – 75° and from +75° to +90°. We calculated the average *c̄*(*lat,t*) and the standard deviation *σ _{c}*(

*lat,t*) of

*c*(

*t*) for all the measurements of each latitude band and of each month: ”

*lat*” is the index that refers to the latitude band and runs from 1 to 17 and ”

*t*” is the index that refers to the month and runs from 1 to 60 (for a total of 5 years). Furthermore, in order to obtain a global value for the stratospheric water vapor trend, we also calculated the average and the standard deviation of

*c*(

*t*) considering for each month all the filtered measurements regardless of latitude. For some months or latitude bands there are no values for

*c̄*(

*lat,t*). This situation occurs in periods for which either nominal MIPAS measurements are not available or the measurements are filtered out due to large errors (see subsection 3.2)

The values of the standard deviation *σ _{c}*(

*lat,t*) are due to several causes of spread: the noise errors that affect the observations, the random term of the incompleteness errors, some instrumental random errors that we discuss in subsection 3.4, the geographic variability due to atmospheric dynamics, the time variability in the month and the dependence on latitude. The error components that are not included in

*σ*(

_{c}*lat,t*) are some systematic instrumental errors that we discuss in subsection 3.4 and the bias term of the incompleteness errors. Therefore, we associate with each value of

*c̄*(

*lat,t*) an error given by the square root of the sum of the variance of the mean plus the squared bias term of the incompleteness error (calculated as the mean of the incompleteness errors of all the measurements in the corresponding latitude band and month).

Figures 3 and 4 show *σ _{c}*(

*lat,t*) as a function of time for each latitude band as well as for all the globe. We tested that a further division of the latitude bands in 5°steps does not decrease the values of

*σ*(

_{c}*lat,t*), showing that the observed spread is dominated by time and longitudinal variability rather than dependence on latitude. From Figs. 3 and 4 we can see that the largest spread of the water vapor concentrations occurs in winter for latitudes higher than 45°in the northern hemisphere and higher than 35°in the southern hemisphere. Very small values of the standard deviation of the order of 0.2 ppmv are obtained when little time and longitudinal variability is present.

In order to estimate the trends of the water vapor, we make the hypothesis that the dependence of *c̄*(*lat,t*) on time can be modeled as the sum of three terms:

*k*(

*lat*) is a constant depending on the latitude band (that characterizes the water vapor content at the beginning of the analyzed time period),

*s*(

*lat,t*) is a periodic function of time with period of one year depending on the latitude band (it describes the seasonal variability) and

*a*(

*lat*) is the slope of the trend that we aim to determine. As we will see below Eq. (19) does not completely represent the time variability of

*c̄*(

*lat,t*), since other atmospheric phenomena with time period different from one year can be present. In general interannual variability, quasi biennial oscillations [29] and processes with periods of several years also affect the observed changes, but the linear component of the trend is the main effect that can be observed in the five year time period and the error with which this trend is fitted will confirm this assumption.

In order to isolate the trend term *a*(*lat*)*t* in Eq. (19) we evaluated the average of *c̄*(*lat,t*) for each month of the year over the five years obtaining a function *c̿*(*lat,m*) (*m* = 1,...,12) which depends on *m* in the same way as *c̄*(*lat,t*) depends on *t* plus a bias term due to the trend. We extended periodically the function *c̿*(*lat,m*) over the five years and then subtracted this function from *c̄*(*lat,t*) obtaining, according to our model, a step function with step width equal to 12 and step altitude equal to 12*a*(*lat*). A fit of this difference with a straight line provided a first estimation of the slope *a*(*lat*) that we used to remove the trend and the bias from *c̿*(*lat,m*). This procedure provided an estimation of the constant part plus the seasonal variability of water vapor content. We subtracted this estimation from *c̄*(*lat,t*) obtaining the final measurement of the trend *a*(*lat*)*t* that was fitted with a straight line. This fit provided a more accurate value of the slope *a*(*lat*). The values obtained for *a*(*lat*) at this second fit are about 8.5% larger than those obtained at the first fit.

In Figs. 5 and 6 we report the estimation of the constant part plus the seasonal variability of water vapor content as a function of time for each latitude band and for all the globe. These figures show that the water vapor VMR is larger at the poles with respect to the equator and that the seasonal variability changes significantly with latitude.

In Figs. 7 and 8 we report the estimation of the trend *a*(*lat*)*t* in percentage (that is divided by the mean of all *c̄*(*lat,t*) over all the months and multiplied by 100) with its best linear fit as a function of time for each latitude band and for all the globe. We can see that the data are well fitted by a straight line in all the latitude bands except in the tropical regions where a strong signature of the quasi biennial oscillation [29] is visible.

In Fig. 9 we report the fitted annual trend 12*a*(*lat*) (black solid line) in percentage (that is divided by the mean of all *c̄*(*lat,t*) over all the months and multiplied by 100) as a function of latitude. The error bars are calculated multiplying the errors obtained from the fit by the square root of the reduced chi-square, that is the chi-square of the fit divided by the number of points minus two (we retrieve two parameters from the fit: the slope and the intercept of the straight line). The fitted annual trend considering all the globe is (0.85±0.08)%yr^{−1}. Trends are determined with errors that are smaller than the fitted values (this conclusion also applies when the systematic errors are considered, see subsection 3.4) and this result justifies the assumption that the observed interannual variability is dominated by a linear variation (see Eq. (19)). Nevertheless the reduced chi-square is significantly larger than one (with values from 7 to 35) showing that the high precision of the measurements is also detecting other causes of variability.

The trend values obtained in the tropical regions are affected by larger errors (as shown by the error bars in Fig. 9) that are most probably due to the presence of the quasi biennial oscillation not included in the model. Figure 9 shows that there is a dependence of the water vapor trend on latitude with values in the northern hemisphere larger than in the southern hemisphere. Neglecting the tropical regions, the mean value of water vapor trend for latitudes higher than 30° is equal to (1.17 ± 0.05)%yr^{−1} in the northern hemisphere and to (0.77 ± 0.035)%yr^{−1} in the southern hemisphere. Therefore, the difference between the percentage values in the two hemispheres is equal to (0.40±0.06)%yr^{−1}. A very smooth dependence on latitude is observed. The Fourier transform of this distribution provides an estimation of the white noise of these points that is smaller than the errors shown by the error bars and estimated on the basis of the geographical variability. This result confirms the high quality of the trend measurements, whose errors tend to be dominated by systematic errors (see subsection 3.4) rather than precision errors.

Since the trend value depends on latitude we make a more accurate estimation of the global trend weighting the latitude dependent values with the fraction of global solid angle seen by each latitude band:

*φ*and

_{max}*φ*are the latitudes delimiting the band. This calculation provides the value of (0.91 ± 0.05)%yr

_{min}^{−1}for the global stratospheric water vapor trend.

#### 3.4. Analysis of the instrumental errors

It is important to evaluate the instrumental calibration errors that can drift during the mission and affect the values of the estimated trends. We analyze the three possible instrumental calibration errors: the frequency calibration error, the error on the instrumental line shape (ILS) and the gain calibration error.

The frequency calibration error affects the spectral axis of the limb spectra measured by MIPAS and can be due to spectral shifts of the diode laser used to sample the interferograms. In order to correct this error a spectral calibration is performed every four elevation scans by comparing measured atmospheric spectra in selected spectral windows with theoretical spectra. The spectral positions of pre-selected lines are compared between the observed spectra and theoretical spectra. The measured spectra are then corrected so that the spectral positions of the compared lines match those of the reference lines. Because of this frequent spectral calibration the frequency calibration error is always very small and appears as a random error in the estimation of the water vapor content and, therefore, it has already been considered in the error estimation presented in subsection 3.3.

In the fit of the measured radiances with the simulated ones the forward model performs a convolution of the simulated spectra with a modeled ILS that should reproduce the real ILS. The modeled ILS is obtained from a model of the instrument and it depends on some parameters whose values were deduced from a fit performed at the beginning of the mission. A change of the real ILS during the mission due to the aging of the instrument (for example a drift of the width of the ILS) can generate an error on the estimated trends. We performed some sensitivity tests changing the values of the parameters characterizing the modeled ILS and estimating the effects of these changes on the retrieved stratospheric water vapor content. From these tests we deduce that the error on the estimated percentage annual water vapor trends due to drifting of the ILS is most probably negligible and in any case less than 0.1%yr^{−1}.

The MIPAS spectra can be affected by a radiometric calibration error meanly due to the effect of the ice contamination on the gain function. In order to eliminate this error a radiometric calibration is performed once a day by observation of an internal black body and the retrieved gain values are updated in the auxiliary data files once a week unless the variation of the gain is greater than 1%. Because of this frequent radiometric calibration the gain error appears as a random error on the retrieved water vapor content and, therefore, it is already considered in the error estimation presented in subsection 3.3. However, the detectors which are used in the long wavelength channels (685 - 1500 cm^{−1}) show a substantial nonlinearity, i.e., their responsivity is dependent on the incident photon flux. This nonlinearity needs to be corrected before the radiometric calibration is applied. The detector response as a function of the incident photon flux has been characterized and parameterized before the ENVISAT launch and the parameters of this characterization are used to perform the nonlinearity correction. A drift of the detector nonlinearity during the mission can induce a drift of the gain error that can affect the estimated stratospheric water vapor trends.

In-flight measurements of the parameters of the nonlinearity correction have allowed to estimate a value of the drift during all the mission of the gain error due to nonlinearity of the detectors of about 2% in band A (that ranges from 685 to 970 cm^{−1}), 1% in band AB (that ranges from 1020 to 1170 cm^{−1}) and 2.5% in band B (that ranges from 1215 to 1500 cm^{−1}), with a sign that determines a compression of the spectrum amplitude as a function of time [30]. We performed a sensitivity test calculating the percentage variation of the retrieved average stratospheric water vapor content for June and December 2009 induced by the estimated drift of the gain error. The results of this sensitivity test (scaled to provide the effect in one year) are shown in Fig. 9 with the blue and red solid lines, respectively. The global variations calculated weighting the latitude dependent values according to Eq. (20) are 0.75%yr^{−1} for June and 0.77%yr^{−1} for December.

These results indicate that the trend values estimated in subsection 3.3 are affected by a bias of about 0.76%yr^{−1}, that, as the blue and red solid lines in Fig. 9 show, is nearly independent of latitude. Combining an uncertainty of 15% on the drift of the gain error with the error due to drifting of the ILS we estimate an error associated to this bias of 0.15%yr^{−1}. Therefore, we conclude that the trend values estimated in subsection 3.3 are affected by a bias of (0.76 ± 0.15)%yr^{−1} which is indicated by the gray area in Fig. 9. Subtracting this bias from the trend values estimated in subsection 3.3 we find (0.15 ± 0.16)%yr^{−1} for the global value, (0.41 ± 0.16)%yr^{−1} for the northern hemisphere and (0.01 ± 0.15)%yr^{−1} for the southern hemisphere.

The considerations made in this subsection show that the drift of the gain error due to nonlinearity of the detectors is the most important source of error that presently drastically reduces the potentiality of the MIPAS measurements for an accurate estimation of stratospheric water vapor trends.

## 4. Conclusion

We developed a new method that uses the MSS to estimate trends of atmospheric parameters from remote sensing measurements. With respect to traditional methods using the retrieved profiles, the proposed method exploits more comprehensively the information included in the observations, avoids the artifacts introduced by the a priori information as well as the loss of information due to interpolation to different vertical grids. We used this new method to estimate the trend of stratospheric water vapor as a function of latitude in the five years from 2005 to 2009 using the MIPAS measurements.

The method provided trend measurements that are very precise (∼5% error), also considering the short time duration of the analyzed period. The estimated trends, however, are affected by a relatively large systematic error induced by the radiometric calibration error of the instrument. The analysis showed that the stratospheric water vapor trend has a dependence on latitude with a positive value of (0.41± 0.16)%yr^{−1} in the northern hemisphere, and less than 0.16%yr^{−1} in the southern hemisphere.

The trend values obtained in this paper show an increase of stratospheric water vapor smaller than that occurred in the period 1980-2000 and in agreement with the small positive trends observed up to 2008 in [14] after the minimum occurred between 2004 and 2006. Our results indicate that the increase is now mainly localized in the northern hemisphere.

The climatological information on the water vapor trend is unfortunately limited by the relatively large radiometric calibration error and by the limited time interval covered by the analyzed data set, but the quality of the results presented in this paper, while providing a good illustrative application of the proposed approach, allows the determination of statistically significant variations in time and as a function of latitude. The good performance of the proposed technique suggests that any atmospheric remote sensing measurement involving an inversion would benefit from a reanalysis with this method.

As a final consideration we stress the importance for a better characterization of the MIPAS radiometric calibration error. Indeed, our analysis revealed that very precise estimates of the trends of the atmospheric parameters can be derived from the MIPAS measurements, but their accuracy is presently drastically limited by the systematic errors.

## Acknowledgments

The study was conducted as part of the CTOTUS project that is supported by the
POR programme of the Tuscany Region within the CReO/FESR objectives
2007-2013
. The forward and inverse models used in this work have been optimized and tested in the frame of the project *Support to MIPAS Level 2 product validation* under
ESA-ESRIN contract N.
21719/08/I-OL.

## References and links

**1. **J. T. Kiehl and K. E. Trenberth, “Earths annual global mean energy budget,” Bull. Am. Meteorol. Soc. **78**, 197–208 (1997). [CrossRef]

**2. **I.M. Held and B.J. Soden, “Water vapor feedback and global warming,” Annu. Rev. Energy Environ. **25**, 441–475 (2000). [CrossRef]

**3. **Intergovernmental Panel on Climate Change (IPCC), *Climate Change 2007: The Physical Science Basis, Contribution of Working Group I to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change*, edited by: S. Solomon, D. Qin, M. Manning, Z. Chen, M. Marquis, K.B. Averyt, M. Tignor, H.L. Miller, (Cambridge University Press, 2007). [PubMed]

**4. **S. Salomon, K. H. Rosenlof, R. W. Portmann, J. S. Daniel, S. M. Davis, T. J. Sanford, and G. K. Plattner, “Contributions of stratospheric water vapor to decadal changes in the rate of global warming,” Science , **327**, 1219–1223 (2010). [CrossRef]

**5. **D. Kley, J. M. Russell, and C. Phillips, “SPARC Assessment of Upper Tropospheric and Stratospheric Water Vapour,” WCRP Report No. 113, WMO/TD Report No. 1043, World Climate Research Programme, Geneva, 325 pp. (2000).

**6. **K. H. Rosenlof, S. J. Oltmans, D. Kley, J. M. Russell III, E. W. Chiou, W. P. Chu, D. G. Johnson, K. K. Kelly, H. A. Michelsen, G. E. Nedoluha, E. E. Remsberg, G. C. Toon, and M. P. McCormick, “Stratospheric water vapor increases over the past half-century,” Geophys. Res. Lett. **28**, 1195–1198 (2001). [CrossRef]

**7. **M. Scherer, H. Vomel, S. Fueglistaler, S. J. Oltmans, and J. Staehelin, “Trends and variability of midlatitude stratospheric water vapour deduced from the re-evaluated Boulder balloon series and HALOE,” Atmos. Chem. Phys. **8**, 1391–1402 (2008). [CrossRef]

**8. **J. M. Russell III, L. L. Gordley, J. H. Park, S. R. Drayson, W. D. Hesketh, R. J. Cicerone, A. F. Tuck, J. E. Frederick, J. E. Harries, and P. J. Crutzen, “The halogen occultation experiment,” J. Geophys. Res. **98**, 10777–10797 (1993). [CrossRef]

**9. **S. J. Oltmans, H. Vomel, D. J. Hofmann, K. H. Rosenlof, and D. Kley, “The increase in stratospheric water vapor from balloon borne, frostpoint hygrometer measurements at Washington, DC, and Boulder, Colorado,” Geophys. Res. Lett. **27**, 3453–3456 (2000). [CrossRef]

**10. **D. J. Seidel, R. J. Ross, J. K. Angell, and G. C. Reid, “Climatological characteristics of the tropical tropopause as revealed by radiosondes,” J. Geophys. Res. **106**, 7857–7878 (2001). [CrossRef]

**11. **X. L. Zhou, M. A. Geller, and M. H. Zhang, “Cooling trend of the tropical cold point tropopause temperatures and its implications,” J. Geophys. Res. **106**, 1511–1522 (2001). [CrossRef]

**12. **G. E. Nedoluha, R. M. Bevilacqua, R. M. Gomez, B. C. Hicks, J. M. Russell III, and B. J. Connor, “An evaluation of trends in middle atmospheric water vapor as measured by HALOE, WVMS, and POAM,” J. Geophys. Res. **108**, 4391–4400 (2003). [CrossRef]

**13. **W. J. Randel, F. Wu, S. J. Oltmans, K. Rosenlof, and G. E. Nedoluha, “Interannual changes of stratospheric water vapor and correlations with tropical tropopause temperatures,” J. Atmos. Sci. **61**, 2133–2148 (2004). [CrossRef]

**14. **A. Jones, J. Urban, D. P. Murtagh, P. Eriksson, S. Brohede, C. Haley, D. Degenstein, A. Bourassa, C. von Savigny, T. Sonkaew, A. Rozanov, H. Bovensmann, and J. Burrows, “Evolution of stratospheric ozone and water vapor time series studied with satellite measurements,” Atmos. Chem. Phys. **9**, 6055–6075 (2009). [CrossRef]

**15. **H. Fischer, M. Birk, C. Blom, B. Carli, M. Carlotti, Tv Clarmann, L. Delbouille, A. Dudhia, D. Ehhalt, M. Endemann, J. M Flaud, R. Gessner, A. Kleinert, R. Koopman, J. Langen, M. Lpez-Puertas, P. Mosner, H. Nett, H. Oelhaf, G. Perron, J. Remedios, M. Ridolfi, G. Stiller, and R. Zander, “MIPAS: an instrument for atmospheric and climate research,” Atmos. Chem. Phys. **8**, 2151–2188 (2008). [CrossRef]

**16. **C. D. Rodgers,
*Inverse Methods for Atmospheric Sounding: Theory and Practice,* Vol. 2 of *Series on Atmospheric, Oceanic and Planetary Physics*
(World Scientific, 2000). [CrossRef]

**17. **B. Carli, P. Raspollini, M. Ridolfi, and B. M. Dinelli, “Discrete representation and resampling in limb-sounding measurements,” Appl. Opt. **40**, 1261–1268 (2001). [CrossRef]

**18. **S. Ceccherini, P. Raspollini, and B. Carli, “Optimal use of the information provided by indirect measurements of atmospheric vertical profiles,” Opt. Express **17**, 4944–4958 (2009). [CrossRef]

**19. **S. Ceccherini, B. Carli, U. Cortesi, S. Del Bianco, and P. Raspollini, “Retrieval of the vertical column of an atmospheric constituent from data fusion of remote sensing measurements,” J. Quant. Spectrosc. Radiat. **111**, 507–514 (2010). [CrossRef]

**20. **S. Ceccherini, U. Cortesi, S. Del Bianco, P. Raspollini, and B. Carli, “IASI-METOP and MIPAS-ENVISAT data fusion,” Atmos. Chem. Phys. **10**, 4689–4698 (2010). [CrossRef]

**21. **S. Ceccherini, B. Carli, E. Pascale, M. Prosperi, P. Raspollini, and B. M. Dinelli, “Comparison of measurements made with two different instruments of the same atmospheric vertical profile,” Appl. Opt. **42**, 6465–6473 (2003). [CrossRef]

**22. **M. Ridolfi, B. Carli, M. Carlotti, T. v. Clarmann, B. M. Dinelli, A. Dudhia, J. M. Flaud, M. Hpfner, P. E. Morris, P. Raspollini, G. Stiller, and R. J. Wells, “Optimized forward model and retrieval scheme for MIPAS near-real-time data processing,” Appl. Opt. **39**, 1323–1340 (2000). [CrossRef]

**23. **P. Raspollini, C. Belotti, A. Burgess, B. Carli, M. Carlotti, S. Ceccherini, B. M. Dinelli, A. Dudhia, J. M. Flaud, B. Funke, M. Hpfner, M. Lpez-Puertas, V. Payne, C. Piccolo, J. J. Remedios, M. Ridolfi, and R. Spang, “MIPAS level 2 operational analysis,” Atmos. Chem. Phys. **6**, 5605–5630 (2006). [CrossRef]

**24. **S. Ceccherini, “Analytical determination of the regularization parameter in the retrieval of atmospheric vertical profiles,” Opt. Lett. **30**, 2554–2556 (2005). [CrossRef]

**25. **S. Ceccherini, C. Belotti, B. Carli, P. Raspollini, and M. Ridolfi, “Technical Note: Regularization performances with the error consistency method in the case of retrieved atmospheric profiles,” Atmos. Chem. Phys. **7**, 1435–1440 (2007). [CrossRef]

**26. **A. Dudhia, V. L. Jay, and C. D. Rodgers, “Microwindow selection for high-spectral-resolution sounders,” Appl. Opt. **41**, 3665–3673 (2002). [CrossRef]

**27. **European Centre for Medium-Range Weather Forecasts, “91 model level definitions,” http://www.ecmwf.int/products/data/technical/model_levels/model_def_91.html (2007).

**28. **J. J. Remedios, R. J. Leigh, A. M. Waterfall, D. P. Moore, H. Sembhi, I. Parkes, J. Greenhough, M. P. Chipperfield, and D. Hauglustaine, “MIPAS reference atmospheres and comparisons to V4.61/V4.62 MIPAS level 2 geophysical data sets,” Atmos. Chem. Phys. Discuss. **7**, 9973–10017 (2007). [CrossRef]

**29. **M. A. Geller, X. L. Zhou, and M. H. Zhang, “Simulations of the interannual variability stratospheric water vapor,” J. Atmos. Sci. **59**, 1076–1085 (2002). [CrossRef]

**30. **M. Birk German Aerospace Center, Remote Sensing Technology Institute, Experimental Methods, Münchner Straße 20, Oberpfaffenhofen-Wessling 82234 (personal communication, 2011).