Direct inversion of circulation and mixing from tracer measurements - Part 1: Method

Abstract

From a series of zonal mean global stratospheric tracer measurements sampled in altitude vs. latitude, circulation and mixing patterns are inferred by the inverse solution of the continuity equation. As a first step, the continuity equation is written as a tendency equation, which is numerically integrated over time to predict a later atmospheric state, i.e., mixing ratio and air density. The integration is formally performed by the multiplication of the initially measured atmospheric state vector by a linear prediction operator. Further, the derivative of the predicted atmospheric state with respect to the wind vector components and mixing coefficients is used to find the most likely wind vector components and mixing coefficients which minimize the residual between the predicted atmospheric state and the later measurement of the atmospheric state. Unless multiple tracers are used, this inversion problem is under-determined, and dispersive behavior of the prediction further destabilizes the inversion. Both these problems are addressed by regularization. For this purpose, a first-order smoothness constraint has been chosen. The usefulness of this method is demonstrated by application to various tracer measurements recorded with the Michelson Interferometer for Passive Atmospheric Sounding (MIPAS). This method aims at a diagnosis of the Brewer-Dobson circulation without involving the concept of the mean age of stratospheric air, and related problems like the stratospheric tape recorder, or intrusions of mesospheric air into the stratosphere.

Atmos. Chem. Phys., 16, 1456314584, 2016 www.atmos-chem-phys.net/16/14563/2016/ doi:10.5194/acp-16-14563-2016 Author(s) 2016. CC Attribution 3.0 License.

Thomas von Clarmann and Udo Grabowski

Karlsruhe Institute of Technology, Institute of Meteorology and Climate Research, Karlsruhe, Germany

Correspondence to: Thomas von Clarmann ([email protected])

Received: 13 April 2016 Published in Atmos. Chem. Phys. Discuss.: 6 June 2016 Revised: 7 October 2016 Accepted: 31 October 2016 Published: 23 November 2016

Abstract. From a series of zonal mean global stratospheric tracer measurements sampled in altitude vs. latitude, circulation and mixing patterns are inferred by the inverse solution of the continuity equation. As a rst step, the continuity equation is written as a tendency equation, which is numerically integrated over time to predict a later atmospheric state,i.e., mixing ratio and air density. The integration is formally performed by the multiplication of the initially measured atmospheric state vector by a linear prediction operator. Further, the derivative of the predicted atmospheric state with respect to the wind vector components and mixing coefcients is used to nd the most likely wind vector components and mixing coefcients which minimize the residual between the predicted atmospheric state and the later measurement of the atmospheric state. Unless multiple tracers are used, this inversion problem is under-determined, and dispersive behavior of the prediction further destabilizes the inversion. Both these problems are addressed by regularization. For this purpose, a rst-order smoothness constraint has been chosen. The usefulness of this method is demonstrated by application to various tracer measurements recorded with the Michelson Interferometer for Passive Atmospheric Sounding (MIPAS). This method aims at a diagnosis of the BrewerDobson circulation without involving the concept of the mean age of stratospheric air, and related problems like the stratospheric tape recorder, or intrusions of mesospheric air into the stratosphere.

Direct inversion of circulation and mixing from tracer measurements Part 1: Method

1 Introduction

In the context of climate change, possible changes in the intensity of the BrewerDobson circulation have become an important research topic. Climate models predict an intensication of the BrewerDobson circulation (Butchart et al., 2006). Engel et al. (2009), however, found a weakly signi-cant slow increase in the mean age of stratospheric air. The latter is dened as the mean time lag between the date of the transition of tropospheric air into the stratosphere and the date when the mixing ratio of a monotonically growing tracer was measured in the air volume under investigation, and its increase suggests a deceleration of the BrewerDobson circulation. These measurements have been challenged as not being representative (Garcia et al., 2011), and global mean age of air measurements by Stiller et al. (2012) suggest that the true picture is not that one-dimensional. Instead, stratospheric age trends vary with altitude and with latitude. The determination of the age of air and its use as a diagnostic of the intensity of the BrewerDobson circulation, however, has its own limitations. First, due to mixing processes, the age of a stratospheric air volume is not unique but characterized by an age spectrum, which has to be considered since the tropospheric growth of SF6 mixing ratios is not strictly linear; some ad hoc assumptions on this spectrum have to be made (Waugh and Hall, 2002). These include the adequacy of the Wald (inverse Gaussian) function for the representation of the age spectrum and the choice of its width parameter. Second, the most suited age tracer, SF6, which has signicant and monotonic growth rates in the troposphere, is not fully inert: it has a mesospheric sink (Hall and Waugh, 1998; Reddmann et al., 2001) and introduces some age uncertainty when mesospheric air subsides into the stratosphere in the polar winter vortex (Stiller et al., 2008). Third, the de-

Published by Copernicus Publications on behalf of the European Geosciences Union.

14564 T. von Clarmann and U. Grabowski: Circulation and mixing from tracer measurements

termination of the mean age relies upon a reference air mass where the age, by denition, is 0. When the age of air concept was introduced, the reference was simply the troposphere, which is well mixed and thus avoids any related complication (Solomon, 1990; Schmidt and Khedim, 1991). Since the age of air has become a model diagnostic, parts of the community have established the stratospheric entry point as a reference (Hall and Plumb, 1994), which makes a difference due to the slow ascent of air through the tropical tropopause layer (Fueglistaler et al., 2009). For model validation, however, this redened age of air is of limited use because no long measured time series of tracer mixing ratios are available there.

Facing these difculties, it is desirable to infer the atmospheric circulation directly from tracer measurements, without going back to the age of air concept. Multiple approaches have been developed to infer wind elds from measured atmospheric state variables. Sequential data assimilation and, in its optimal form, the extended Kalman lter approach (e.g., Ghil and Malanotte-Rizzoli, 1991; Ghil, 1997), calculates the optimal average of the forecasted meteorological variables for the time of the observation and the observed meteorological variables themselves and uses this average to initialize the next forecast step. The wind eld is calculated by a dynamical model. This method involves the generalized inversion of the observation operator where the forecast is used as a constraint. In contrast, so-called1 variational data assimilation minimizes the residual between the forecasted and the measured atmospheric state variables by optimally adjusting the initialization of the forecast via inversion of an adjoint forecast model, constrained by some background state (Thompson, 1969). Both approaches rely on dynamical models2 and are suited to infer the most probable atmospheric state variables rather than the wind eld, which is a by-product of the assimilation. The wind eld or atmospheric circulation can also be inferred directly by kinematic methods from tracer measurements. Such methods rely solely on the continuity equation, do not involve a dynamic model, and thus do not depend on any ad hoc parameterization of effects which are either not resolved by the discrete model, computationally too expensive for explicit modeling, or simply not well understood. While this work is targeted primarily at an assessment of the BrewerDobson circulation, its applicability is much wider and includes stratospherictropospheric exchange, the mesospheric overturning circulation, and other areas. Early approaches to infer the circulation from tracer measurements include Holton and Choi (1988) as well as Salby and Juckes (1994) who used approaches which share several ideas with ours. The direct inversion of wind speeds

1The term so-called is used here because it is challenged that this method is really variational in the context of discrete variables (Wunsch, 1996, p. 368).

2This statement refers to meteorological data assimilation. Chemical data assimilation uses chemistry transport models.

from tracer measurements in a volcano plume has, e.g., been suggested by Krueger et al. (2013), however without consideration of mixing. The continuity equation including diffusion terms has been exploited by Wofsy et al. (1994) for the assessment of diffusion of stratospheric aircraft exhaust.

In this paper we present a method to infer two-dimensional (latitudealtitude) circulation and mixing coefcients from subsequent measurements of inert tracers. The application of this method, i.e., the inference of the BrewerDobson circulation from global SF6 distributions (Stiller et al., 2008, 2012) measured with the Michelson Interferometer for Passive Atmospheric Sounding (MIPAS), is presented in a companion paper. In order to avoid the reader not seeing the forest for the trees, we give a short overview of our method in Sect. 2. The prediction of pressure and tracer mixing ratio elds on the basis of the continuity equation and related error estimation is described in Sect. 3. The estimation of circulation and mixing coefcients by the inversion of the continuity equation is presented in Sect. 4. In Sect. 5, the applicability of the method and the need for further renements is discussed critically. The benets of the method are discussed in Sect. 6. The paper concludes with recommendations on how these results should be used and with an outlook for future work (Sect. 7). Changes in the BrewerDobson circulation during 20022012, i.e., the MIPAS mission, are currently investigated by means of this method and will be published in a subsequent paper.

2 General concept

Knowing the initial state of the atmosphere in terms of mixing ratio and air density distributions, wind speed, and mixing coefcients at each grid point, a future atmospheric state can be predicted with respect to the distribution of any inert tracer. This procedure we call the forward problem. If no ideal tracers are available, source and sink terms of related species have to be included in the forward model. The goal of this work is to invert the forward model in order to infer the circulation and mixing coefcients from tracer measurements by the minimization of the residual between the predicted and measured atmospheric state. This approach is complementary to free-running climate models because it makes no assumptions about atmospheric dynamics except for the validity of the continuity equation. It is further considered more robust than age-of-air analysis (Stiller et al., 2012) because it does not depend on a reference point where the age is assumed 0, nor does it require knowledge of the history of an air parcel.

Our concept involves the following operations. First, a general solution of the forward problem is formulated (Sect. 3). The forward problem is the solution of the prediction equation as a function of the initial atmospheric state for given winds and mixing coefcients. For our chosen solver, which involves the MacCormack (1969) integration scheme

Atmos. Chem. Phys., 16, 1456314584, 2016 www.atmos-chem-phys.net/16/14563/2016/

T. von Clarmann and U. Grabowski: Circulation and mixing from tracer measurements 14565

for the solution of the transport problem (Eqs. 510), the relevant dependencies of the nal state on the initial state are reported in Sect. 3.2 (Eqs. 1526). The formulation in matrix notation (Eqs. 2728) allows the easy treatment of multiple successive time steps (Eq. 29) and an easy estimation of the prediction error via generalized Gaussian error estimation (Eq. 30). As a next step, the dependence of the predicted state on winds and mixing coefcients is estimated for a given initial state. This is achieved by the differentiation of the solution of the prediction equation with respect to winds and mixing coefcients (Eqs. 3476). These partial derivatives form the Jacobian matrix of the problem, with which the estimation of winds and mixing coefcients can be reduced to a constrained least squares optimization problem where the inversely variance-weighted residual between the predicted atmospheric state and the respective measured atmospheric state is minimized. The latter step involves the generalized inverse of the Jacobian matrix (Eqs. 7890).

3 The forward problem

The forward model reads the measured atmospheric state at time t and predicts the atmospheric state (number density of air, c, and volume mixing ratios, vmr) at time t +[Delta1]t for given

wind vectors and mixing coefcients representing the time interval [notdef]t;t + [Delta1]t[notdef] by solving the continuity equation. The

continuity equation allows us to calculate the local tendencies of the number densities and volume mixing ratios. These local tendencies @@t and @vmr@t are then integrated over time to give the new number densities and mixing ratios.

3.1 The continuity equation

The local change in number density of air is in spherical coordinates (for all auxiliary calculations, see Supplement 1):

@

@t =

1 r

@vmrg

@vmrg

@vmrg

@ (2)

w @vmrg @z +

@t =

Sg

u

r cos

@

v r

1 r2

@ @

[bracketleftbigg]

K cos2

@vmrg @

[bracketrightbigg]

+

1r2 cos

@ @

K cos @vmrg @

[bracketrightbigg]

+

1 r2

@ @z

r2Kz @vmrg @z

,

where vmrg is the volume mixing ratio of species g, K is the zonal diffusion coefcient, K is the meridional diffusion coefcient, Kz is the vertical diffusion coefcient, and Sg is the production/loss rate of species g in terms of number density over time (e.g., Brasseur and Solomon, 2005, Eq. 3.46b, and Jones et al., 2007).

Since we are only interested in a two-dimensional representation of the atmosphere in altitude and latitude coordinates, zonal advection and mixing terms are ignored in Eqs. (1)(2). In this two-dimensional representation, all atmospheric state variables represent zonal mean values. The kinematic variables, viz. the velocities and mixing coefcients, have to be reinterpreted because they do not represent merely the zonally averaged velocities and mixing coefcients. Instead, they include also eddy transport and diffusion terms, and their interpretation is less unique than one might hope because it depends on the denition of the kinematic variables and the approximations used (see Appendix A for details of the interpretation of the kinematic quantities). The local change in number density of air in a two-dimensional atmosphere thus is

@

@t =

1 r

@v

@ +

v

r tan()

2w

r , (3)

and the local change in vmrg is calculated as

@vmrg

@t =

Sg

@w

@z

@v

@ +

v

r tan()

@w

@z

2w

r (1)

1r cos()

@u

@ ,

v r

@vmrg

@ w

@vmrg

@z (4)

+

1r2 cos

@ @

K cos @vmrg @

[bracketrightbigg]

+

where t is time, is longitude, is latitude, z is altitude above surface, r = rE + z, rE is the radius of Earth,

u = (rE + z)cosd /dt, v = (rE + z)d/dt, and w = dz/dt.

Here the shallowness approximation (Hinkelmann, 1951; Phillips, 1966, quoted after Kasahara, 1977), which simplies the equations using the assumption that z is much smaller than rE and which is, often implicitly, used in the usual textbooks on atmospheric sciences (e.g., Brasseur and Solomon, 2005, their Eq. 3.46a), is intentionally not used for reasons which will become clear in Sect. 3.2.

The local change in the volume mixing ratio of gas g can be calculated from known velocities and mixing coefcients as well as source/sink terms as

1 r2

@ @z

r2Kz @vmrg @z

.

3.2 Integration of tendencies

The integration of Eqs. (3)(4) is performed numerically for time steps of [Delta1]tp. For practical reasons, processes (advection, diffusion, sinks) are split; i.e., the tendencies due to these three classes of processes are integrated independently. The time steps [Delta1]tp used for the integration are chosen to be smaller than the time difference [Delta1]t between two measurements, in order not to clash with the Courant limit

www.atmos-chem-phys.net/16/14563/2016/ Atmos. Chem. Phys., 16, 1456314584, 2016

14566 T. von Clarmann and U. Grabowski: Circulation and mixing from tracer measurements

(Courant et al., 1952). In the following, we call [Delta1]tp a micro time increment and [Delta1]t macro time increment. The atmospheric state after a macro time increment is predicted by successive prediction over the micro time increment. In the following, index i designates time t, i + 1 designates the

time t + [Delta1]tp, etc., and I designates the time after the nal

micro time increment, i.e., the next macro time increment.

For the discrete integration of the advection part of the tendencies, the MacCormack (1969) method is used in a generalized multidimensional version similar to the one described by (Perrin and Hu, 2006). This is a predictorcorrector method. For a general state variable c(t,x,y) = ci(x,y) at

location (x,y), and time t with e(c) and f (c) being functions of c, an equation of the form

@c

@t +

@e(c)

@x +

@f (c)

and the corrected prediction for then gives

i+1(,z) =

1

2r2 cos() [notdef] [bracketleftbigg]

i(,z)r2 cos() (10)

+

hi+1(,z)r2 cos()

i

[Delta1]tp

[Delta1]

[bracketleftbig][bracketleftbig]

i+1(,z)rv(,z)cos()

i+1( [Delta1],z)rv( [Delta1],z)cos( [Delta1])

i

[bracketleftBig][bracketleftBig]

i+1(,z)r2w(,z)cos()

[bracketrightbig]

[Delta1]tp

[Delta1]z

@y = 0 (5)

is solved by preliminary predictions of the state variable as a rst step: x is

c i+1(x,y) = ci(x,y) (6)

[Delta1]tp

[Delta1]x (ei(x + [Delta1]x,y) ei(x,y))

[Delta1]tp

hi+1(,z [Delta1]z)(r [Delta1]z)2w(,z [Delta1]z)cos()[bracketrightBig] [bracketrightBig][bracketrightbigg] .

For the local change in mixing ratio, operator splitting is performed. The horizontal and vertical advective parts of the continuity equation for mixing ratios in two dimensions are transformed into the following MacCormack-integrable forms:

@ r vmrv @t

adv.horiz=

@vmrg

@ (11)

and

"@ vmrgw @t

[bracketrightBigg]adv.vert =

@vmrg

@z , (12)

respectively.

For the diffusive component we use simple Eulerian integration:

vmrg;i+1(,z) vmrg;i(,z)

diff (13)

=

[Delta1]y (fi(x,y + [Delta1]y) fi(x,y)).

These are then used in a subsequent correction step which gives the nal prediction

ci+1(x,y) =

1

2

hci(x,y) + c i+1(x,y) (7)

[Delta1]tp

[Delta1]x e(c i+1,x,y) e(c i+1,x [Delta1]x,y)

[parenrightbig]

[Delta1]tp

[Delta1]y f (c i+1,x,y) f (c i+1,x,y [Delta1]y)

[parenrightbig][bracketrightBig]

.

[Delta1]tp

2r2([Delta1])2 cos()

[notdef]

[bracketleftbigg][parenleftBigg]K

(,z) + K( + [Delta1])

[parenrightbig]

Application to the continuity equation in spherical coordinates requires the reformulation of Eq. (3) (see, e.g., Chang and St.-Maurice, 1991):

@r2 cos()@t =

@r v cos()@

cos( +

[Delta1]

2 )

(vmrg;i( + [Delta1],z) vmrg;i(,z))

K(,z) + K( [Delta1])

[parenrightbig]

cos(

[Delta1]

2 )

@r2 w cos()@z . (8)

The predictor of r2 cos() is then calculated as

[r2i+1(,z))cos()[notdef] = r2i(,z)cos() (9)

[Delta1]tp

[Delta1] (ri( + [Delta1],z)cos( + [Delta1])v( + [Delta1],z)

ri(,z)cos()v(,z))

[Delta1]tp

[Delta1]z

(r + [Delta1]z)2i(,z + [Delta1]z) w,z+[Delta1]z cos()

r2i(,z) w,z cos()[parenrightBig]

,

[notdef] (vmrg;i(,z) vmrg;i( [Delta1],z))

[bracketrightbigg]

(r +

[Delta1]z2 )2

[notdef] (Kz(,z) + Kz(,z + [Delta1]z))

[notdef] (vmrg;i(,z + [Delta1]z) vmrg;i(,z))

r

[Delta1]z 2

2[notdef] (Kz(,z) + Kz(,z [Delta1]z))

[notdef] (vmrg;i(,z) vmrg;i(,z [Delta1]z))

.

+

[Delta1]tp

2r2([Delta1]z)2

Atmos. Chem. Phys., 16, 1456314584, 2016 www.atmos-chem-phys.net/16/14563/2016/

T. von Clarmann and U. Grabowski: Circulation and mixing from tracer measurements 14567

Sinks of the species considered here are treated as uni-molecular processes (see, e.g., Brasseur and Solomon, 2005, their Eq. 2.27d) and integrated as

g;i+1 = g;iekg[Delta1]tp, (14)

where kg is the sink strength of the gas g.

The abundance of gas g after time step [Delta1]tp is then simply the sum of the increments due to horizontal and vertical advection, diffusion, and chemical losses.

Admittedly, there exist more elaborate advection schemes than the one used here. However, the need to provide the Jacobians needed in Sects. 3.34 justies a reasonable amount of simplicity. Further, numerical errors cannot easily accumulate because after each time step [Delta1]t, the system is reinitialized with measured data.

Since we do not have a closed system but have mass exchange and mixing with the atmosphere below the lowermost model altitude and above the uppermost altitude, the atmospheric state is not predicted for the lowermost and uppermost altitudes. Prediction is only possible from the second altitude from the bottom up to the altitude below the uppermost one. Henceforth, we call this restricted altitude range the nominal altitude range. Instead, the atmospheric state of the uppermost and lowermost altitude is estimated by the linear interpolation of measured values at times t and t + [Delta1]t

and used as a boundary condition for prediction within the nominal altitude range. Alternatively, derivatives at the border can be approximated by asymmetric difference quotients.

We use the following convention. Atmospheric state variables are sampled on a regular latitudealtitude grid. For some grid points, no valid measurements may be available but we assume that, for each state variable, we have a contiguous subset of this grid with valid measurements. For state variable g, we have a total of Jg valid measurements within the nominal altitude range, each denoted by index j. A state variable in this context is either air density (g = 0)

or the mixing ratio of one species vmrg. The nominal altitude range at latitude is the altitude range where, for each grid point, a valid measurement is available at the grid point itself and for its northern, southern, upper, lower, and diagonal neighbors. The use of asymmetric difference quotients can be emulated by extrapolation, generating articial border values, which guarantee that each grid point within the nominal altitude and latitude range has all required neighbors. The availability of neighbor values is necessary to allow the calculation of numerical derivatives of the state variable with respect to latitude and altitude. We therefore have Kg border elements of each quantity g, each denoted by index k. This implies, for each state variable g, a total of Lg = Jg + Kg

grid points with indices l.

3.3 Integration in operator notation

For further steps (error propagation and the solution of the inverse problem), it is convenient to rewrite the prediction

of air density and mixing ratios in matrix notation. For this purpose, we differentiate the predicted air densities (Eq. 10) and mixing ratios (Eqs. 1113) with respect to air density and mixing ratios of the gases under assessment at all relevant locations. The sensitivities of the densities of the rst predictive step with respect to the initial densities at the same latitude and altitude are

@i+1(,z)

@i(,z) (15)

=

1

2

2 [Delta1]tp [Delta1]

v(,z) r

[Delta1]tp

[Delta1] [notdef]

v(,z)

r +

[Delta1]tp

[Delta1]z w(,z) [parenrightbigg]

+

[Delta1]tp

[Delta1]

v( [Delta1],z)v(,z)

r2

[bracketrightbigg]

[Delta1]tp

[Delta1]z

w(,z)

[Delta1]tp

[Delta1] [notdef]

v(,z)

r +

[Delta1]tp

[Delta1]z w(,z) [parenrightbigg]

+

[Delta1]tp

[Delta1]z w(,z [Delta1]z)w(,z)

[bracketrightbigg][bracketrightbigg]

.

We further differentiate predicted air densities with respect to air densities at the adjacent southern latitude but the same altitude.

@i+1(,z) @i( [Delta1],z) =

1

2

[Delta1]tp

[Delta1] [notdef]

v( [Delta1],z)

r [notdef]cos( [Delta1]) cos()

(16)

1 +

[Delta1]tp

[Delta1] [notdef]

v( [Delta1],z)

r +

[Delta1]tp

[Delta1]z w( [Delta1],z)

[parenrightbigg][bracketrightbigg]

The derivative of the predicted air densities with respect to air densities at the adjacent northern latitude but the same altitude is

@i+1(,z) @i( + [Delta1],z) =

1

2

[Delta1]tp

[Delta1] [notdef]

v( + [Delta1],z)

r [notdef]cos( + [Delta1]) cos()

(17)

1 +

[Delta1]tp

[Delta1] [notdef]

v(,z)

r +

[Delta1]tp

[Delta1]z w(,z)

[parenrightbigg][bracketrightbigg] .

As a next step we differentiate predicted air densities with respect to the initial air densities at the next higher altitude but the same latitude.

@i+1(,z) @i(,z + [Delta1]z) =

1

2

[Delta1]tp

[Delta1]z [notdef]

(r + [Delta1]z)2r2 w(,z + [Delta1]z) (18)

[notdef]

1 +

[Delta1]tp

[Delta1] [notdef]

v(,z)

r +

[Delta1]t

[Delta1]z w(,z)

[parenrightbigg][bracketrightbigg] .

The derivative of the predicted air densities with respect to the initial air densities at the next lower altitude but the same latitude is

@i+1(,z) @i(,z [Delta1]z) =

1

2

[Delta1]tp

[Delta1]z w(,z [Delta1]z)

(r [Delta1]z)2r2 (19)

[notdef]

1 +

[Delta1]tp

[Delta1]

v(,z [Delta1]z)

r [Delta1]z +

[Delta1]tp

[Delta1]z w(,z [Delta1]z)

[parenrightbigg][bracketrightbigg] .

www.atmos-chem-phys.net/16/14563/2016/ Atmos. Chem. Phys., 16, 1456314584, 2016

14568 T. von Clarmann and U. Grabowski: Circulation and mixing from tracer measurements

Finally we differentiate the predicted air densities with respect to the initial air densities at the adjacent southern latitude and higher altitude,

@i+1(,z) @i( [Delta1],z + [Delta1]z)

(20)

@vmri+1(,z) @vmri( [Delta1],z)

(24)

=

v(,z) 2r [notdef]

[Delta1]tp

[Delta1]

1 +

[Delta1]tp

[Delta1] [notdef]

v( [Delta1],z)

r

[parenrightbigg]

+

=

1

2

v( [Delta1],z)

r [notdef]

[Delta1]tp

[Delta1] [notdef]

[Delta1]tp

[Delta1]z [notdef]

(r + [Delta1]z)2 r2

[Delta1]tp

2r2([Delta1])2 cos()

[notdef]

cos( [Delta1])

[notdef] cos() w( [Delta1],z + [Delta1]z)[bracketrightbigg]

,

K(,z) + K( [Delta1],z)

cos

[Delta1] 2

;

and vice versa

@i+1(,z) @i( + [Delta1],z [Delta1]z)

@vmri+1(,z) @vmri(,z + [Delta1]z)

(25)

(21)

=

1

2 [notdef] w(,z) [notdef]

[Delta1]tp

[Delta1]z

1

[Delta1]tp

[Delta1]z w(,z) [parenrightbigg]

+

=

1

2

w(,z [Delta1]z)

[Delta1]tp

[Delta1]z [notdef]

[Delta1]tp

[Delta1] [notdef]

(r [Delta1]z)2 r2

[notdef]

[Delta1]tp

2r2([Delta1]z)2 (r +

[Delta1]z

2 )2[parenleftbigg]

Kz(,z) + Kz(,z + [Delta1]z)

;

cos( + [Delta1])

cos() [notdef]

v( + [Delta1],z [Delta1]z)

r [Delta1]z

,

@vmri+1(,z) @vmri(,z [Delta1]z)

(26)

where i is the index of the time increment and where [notdef][Delta1]

and z [notdef] [Delta1]z refer to the adjacent model grid points in latitude

and altitude, respectively.For mixing ratios, the respective derivatives are

@vmri+1(,z)

@vmri(,z) (22)

= 1

[Delta1]tp [Delta1]

=

[Delta1]tp

[Delta1]z [notdef]

1

2 [notdef] w(,z)[parenleftbigg]

1 + w(,z [Delta1]z)

[Delta1]tp

[Delta1]z

[parenrightbigg]

+

[Delta1]tp

2r2([Delta1]z)2

r

[Delta1]z 2

2 Kz(,z) + Kz(,z [Delta1]z)

,

where Loss(month,,z) is the relative loss rate in the respective month at latitude and altitude z. These derivatives are simplications in the sense that they do not consider the full chemical Jacobian but assume instead that the source strength depends on no other concentration than the actual concentration of the same species. For the typical long-lived tropospheric source gases considered here, like SF6 or CFCs, this assumption is appropriate. Pre-tabulated loss rates are used which have been calculated by locally integrating loss rates over an entire month at a time resolution adequate to resolve the diurnal cycle. From the monthly losses, the Loss(month,,z) values, which are the contribution of losses to the partial derivatives of the local mixing ratios with respect to the initial local mixing ratios, are calculated as the secant of the local decay curve.

With these expressions, the prediction of air density and volume mixing ratio can be rewritten in matrix notation for a single micro time increment. This notation simplies the estimation of the uncertainties of the predicted atmospheric state and the inversion of the prediction equation. In matrix notation, the prediction reads

i+1 =

2[notdef]

v(,z)

r2 [notdef]1 2

v(,z) + v( [Delta1],z)

[bracketrightbig]

[Delta1]tp [Delta1]z

2[notdef] w(,z) [notdef]1 2

w(,z) + w(,z [Delta1]z)

[bracketrightbig]

[Delta1]tp

2r2([Delta1])2 cos()

[bracketleftbigg][parenleftbigg]

K(,z) + K( + [Delta1],z)

[parenrightbigg]

[notdef] cos

+

[Delta1]2[parenrightbigg]

+

K(,z) + K( [Delta1],z)

cos

[Delta1]2[parenrightbigg][bracketrightbigg]

[Delta1]tp

2r2([Delta1]z)2

[bracketleftbigg][parenleftbigg]

r +

[Delta1]z

2

2 Kz(,z) + Kz(,z + [Delta1]z) [parenrightbigg]

+

r

[Delta1]z 2

2 Kz(,z) + Kz(,z [Delta1]z) [parenrightbigg][bracketrightbigg]

Loss(month,,z)[Delta1]tp; @vmri+1(,z)

@vmri( + [Delta1],z)

(23)

=

[Delta1]tp

[Delta1] [notdef]

v(,z) 2r

1

[Delta1]tp

[Delta1] [notdef]

v(,z)

r

[parenrightbigg]

+

0 @

I;k=1,K0

i+1;k=1,K0

i+1;j=K0+1,L0

1 A

= D;ii

[Delta1]tp

2r2([Delta1])2 cos()

[notdef]

=

IK 0 0

Wi 0 0 D,nom

A,(27)

Atmos. Chem. Phys., 16, 1456314584, 2016 www.atmos-chem-phys.net/16/14563/2016/

0 @

1 A

0 @

I;k=1,K0

i;k=1,K0

i;j=K0+1,L0

1

K(,z) + K( + [Delta1],z)

cos

+

[Delta1] 2

;

T. von Clarmann and U. Grabowski: Circulation and mixing from tracer measurements 14569

where

D;i is the (L0 +K0)[notdef](L0 +K0) Jacobian matrix of air

density for time increment i, i.e., the sensitivities of the prediction with respect to the initial state, @ci+1,m@ci,n (here m and n run over the model grid points);

IK is K0 [notdef] K0 identity matrix; 0 are zero sub-matrices

of the required dimensions;

Wi is a K0 [notdef] 2K0-dimensional interpolation matrix;

D,nom is a J0 [notdef] L0 Jacobian containing the partial

derivatives @i+1;j /@i;l, applied to the nominal alti

tude range;

I;k=1,K0 is the K0-dimensional vector of air densities

in the border region after the nal time step, i.e., for the time of the next measurement;

i;k=1,K0 is the K0-dimensional vector of air densities

in the border region at the current time step as resulting from interpolation in time; and

i;j=K0+1,L0 is the J0-dimensional vector of air densi

ties in the nominal region at the current time step as resulting from integration according to the MacCormack scheme as described above.

The operation of these sub-matrices is illustrated in the uppermost three (violet, green, blue) blocks of Fig. 1.

Since the source term depends on air density, the integration in matrix notation for vmr requires simultaneous treatment of vmr and air density, and we get the following, using notation accordant with air density:

[parenleftbigg]

i+1

vmri+1

[parenrightbigg] =

0

B

B

B

B

B

B

@

I;k=1,K0

i+1;k=1,K0

i+1;j=K0+1,L0

vmrg;I;k=1,Kg vmrg;i+1;k=1,Kg vmrg;i+1;j=Kg+1,Lg

1

C

C

C

C

C

C

A

[parenrightbigg] =[parenleftBigg][parenleftBigg]

1

Yi=I

Di

[parenrightBigg][parenleftbigg]

0 vmr0

. (29)

= Di

[parenleftbigg]

i vmrg;i

[parenrightbigg]

3.4 Prediction errors

Let S0 be the L [notdef] L covariance matrix describing the uncer

tainties of all involved measurements 0 and vmr0, with diagonal elements s0;l,l = 20;l,l and L =

=

0

B

B

@

D;i 0 0 0 0 IK 0 0

0 Wi 0

Dg,nom

1

C

C

A

0

B

B

B

B

B

B

@

I,k=1,K0

i;k=1,K0

i;j=K0+1,L0

vmrg;I,k=1,Kg vmrg;i,k=1,Kg vmrg;i,j=Kg+1,Lg

1

C

C

C

C

C

C

A

,

Pgases0 Lg. We assume that these measurement errors in the state variables used for the prediction are the only relevant error sources. With S0 and

Q1i=I Di available, generalized Gaussian error propagation for I and vmrI can be easily formulated as

SI =

(28)

where Di is the total Jacobian with respect to air densities and all involved gas mixing ratios, and where g runs over all gases. Note that

1. the Jacobian Di is time-dependent because it includes sub-matrices controlling the interpolation between the initial time and the end time. In the case of vmr, a further time dependence is introduced by the time-dependent source function.

2. the rst row of the Jacobian matrix includes identity IK because the prediction is not supposed to change the measured I and vmrI at the end of the macro time increment. This value is used to construct the boundary condition. Row is here written in quotes because the elements of this row are matrices in themselves. The introduction of unity Jacobian elements is necessary because Eqs. (27)(28) are autonomized, originally non-autonomous systems of differential equations.

3. W is used to interpolate the boundary state between the initial time of the micro time interval, t+(i1)[Delta1]tp, and

the time at the end of the time interval t +[Delta1]t to give the

atmospheric state at the border region at time t + i[Delta1]tp.

4. the Jacobian sub-matrices D,nom and Dg,nom are used to predict the atmospheric state in the nominal range after one further micro time increment from the atmospheric state at the current time and the boundary condition. Its elements are described in Eqs. (15)(21) and(22)(26). The part of Dg,nom which refers to the border mixing ratios (vmrg;I,k=1,Kg) is zero. The dimension of

Dg,nom is

PgJg [notdef] (L0 + K0 + Pg(Lg + Kg)).5. no simple mapping mechanism between the eld of atmospheric state variables sampled at latitudes and altitudes and the vectors and vmrg is provided because the elds are irregular in the sense that the number of relevant altitudes is latitude-dependent. Pointer variables have to be used instead.

The matrix structure is exemplied in Fig. 1. For the macro time increment [Delta1]t, we get

[parenleftbigg]

I vmrI

1 !S01

Yi=I

Di

!T. (30)

Even if S0 is diagonal, i.e., the initial errors are assumed to be uncorrelated, error propagation through the forward model will generate non-zero error covariances in SI representing the atmospheric state at time t +[Delta1]t. SI will be needed in the

inversion of circulation and mixing coefcients described in Sect. 4.

www.atmos-chem-phys.net/16/14563/2016/ Atmos. Chem. Phys., 16, 1456314584, 2016

14570 T. von Clarmann and U. Grabowski: Circulation and mixing from tracer measurements

D ;J

L L

L0 g g

r

.. r

i;1

i;K

i;L

r.

.r

I;1

I;K0

0

i;K +1

0

0

g

IK

W

i

r.

.r

r

I;1

i;K

IK

I;K

g

W

vmr

vmr

vmr

vmr

vmr

.

i;1

i;K +1

.

i

vmr

D

g

;J

.

.

g

i;Lg

g

I;1

i;K

vmr

vmr

vmr

vmr

vmr

IK

.

I;K

i;1

i;L

W

.

i

vmr

D

g

i;K +1

g

g

g

;J

.

K0

K J0 Kg Kg J Kg Kg Jg

0

g

Figure 1. Matrix structure of the right-hand part of Eq. (28). Colors indicate which matrix blocks operate on which part of the input vector. An example with two gases is shown.

3.5 A note on nite-resolution measurements

The measurements used are not a perfect image of the true atmospheric state but contain some prior information.In the case of the data provided by the Institute of Meteorology and Climate Research (IMK), a priori proles are usually set to zero, and the constraint is built with a Tikhonov-type rst-order nite differences smoothing constraint (see von Clarmann et al., 2009). That means that, besides the mapping of measurement and parameter errors, the only distortion of the truth via the retrieval is reduced altitude resolution; no other effect of the prior information is to be considered. Usually, any comparison between modeled and measured elds requires application of the averaging kernels of the retrieval to the model data in order to account for the smoothing by the constraint of the retrieval (assuming that the model grid is much ner than the resolution of the retrieval).

In our case, the situation is different. The model is initialized with measurements of reduced altitude resolution, and the elds predicted by the model are then compared to measurements of the same altitude resolution. It is fair to assume that the model does not dramatically change the altitude resolution of the proles, and thus comparable quantities are compared when the residuals between predicted and measured atmospheric state are evaluated.

3.6 A note on numerical mixing

Let the initial mixing ratio eld be homogeneous except one point with delta-type excess mixing ratio. Assume further a homogeneous velocity eld and zero mixing coefcients. If the velocity is such that the position of the excess mixing ratio is displaced during [Delta1]t by a distance which is not equal to an integer multiple of the grid width, then the resulting distribution will no longer be a delta-type distribution but will be smoothed. We refer to the widening of the delta peak as numerical mixing. The MacCormack transport scheme is designed to ght this diffusivity but some higher-order effects may still survive. One might think that, during the inversion, the widening is misinterpreted as mixing, leading to mixing coefcients that are too large.

Again, in our case, the situation is different. The widening does not accumulate over the [Delta1]tp time steps because we rst calculate the operator

Q1i=I Di, which is applied only once to the initial eld, which avoids the accumulation of numerical mixing over time steps. Still one widening process as described above can occur, when the forward model leads to a position of the new peak which cannot be represented in the grid chosen. However, since the gas distributions vmrg;I

at the end of time step [Delta1]t are sampled on the same grid, the maximum in the real atmosphere would be widened in the same way and there would be no residual the inversion would

Atmos. Chem. Phys., 16, 1456314584, 2016 www.atmos-chem-phys.net/16/14563/2016/

T. von Clarmann and U. Grabowski: Circulation and mixing from tracer measurements 14571

try to get rid of by increasing the mixing coefcients. And the next time step [Delta1]t is initialized again with measured data, which also excludes the accumulation of numerical mixing effects.

These considerations aside, there are other numerical artifacts. These are related to the numerical evaluation of partial derivatives of the state variables in our transport scheme chosen. Particularly in the case of delta functions in the state variable eld, these cause side wiggles behind and smearing in front of the transported structure. To keep these artifacts small, it is necessary to set the spatial grid ne enough that every structure is represented by multiple grid points.

4 The inverse problem

For convenience, we combine the variables of the initial atmospheric state and the predicted state at the end of the macro time interval, respectively, into the vectors

ex0 =

[parenleftbigg]

0 vmr0

Eq. (28) by application of the product rule:

ef n =

I

Xi=1

[bracketleftBigg][parenleftBigg]i+1[productdisplay]

[parenrightBigg][parenleftbigg]

@Di

@qn

[parenrightbigg][parenleftBigg]

1

Yk=i1 Dk

[parenrightBigg]

ex0

#, (35)

k=I

Dk

where the tilde symbol in

ef n indicates that the vectors resulting from Eq. (35) still include the border elements which have to be discarded to obtain f n. The quantity

ef nis more efciently computed using the following recursive scheme, where

ef l,i is the respective column of the Jacobian after micro time step i:

ef n,i = Dif n,i1 +

@Di

@qn

1

Yk=i1 Dk

[parenrightBigg]

ex0. (36)

With the argument of D specifying the column of the D matrix such that Dc,i(,z) relates i+1(,z) to i(,z),

D,i( [notdef] [Delta1],z) relates i+1(,z) to i( [notdef] [Delta1],z), and

D,i(,z [notdef] [Delta1]z) relates i+1(,z) to i(,z [notdef] [Delta1]z), and for

vmr accordingly, the entries of Di relevant to v are

@ @i+1(,z)@i(,z)

@v(,z) =

[Delta1]tp

2[Delta1] (37)

[notdef]

[parenrightbigg] =

(

ex0;1...

ex0;L)T (31)

and

exI =

I vmrI

[parenrightbigg] =

(

exI;1...

exI;L)T. (32)

The related subsets of

[Delta1]tp

[Delta1] [notdef]2v(,z) + v( [Delta1],z)

r2 + 2[Delta1]tp

[Delta1]z [notdef] w(,z) r

,

ex0 and

exI which contain only state variables in the nominal altitude range but not those in the border region are x0 = (x0;1...x0;J )T and xI =

(xI;1...xI;J )T, respectively. The reason why the distinction

between

ex and x is made is that, contrary to the prediction step, for the inversion vector elements related to the interpolation of values in the border region are no longer needed. Further, we combine the elds of meridional and vertical wind components and mixing coefcients into the vector

q =

0

B

B

@

@ @i+1(+[Delta1],z)@i(+[Delta1],z)

@v(,z) =

1

2 [notdef] [parenleftbigg]

[Delta1]tp

[Delta1]

2[notdef]

v( + [Delta1],z)

r2 , (38)

@ @i+1(+[Delta1],z)@i(,z)

@v(,z) =

1

2r [notdef]

[Delta1]tp

[Delta1] [notdef]

cos() cos( + [Delta1])

(39)

[notdef]

1 + 2[Delta1]tp

[Delta1] [notdef]

v(,z)

r +

[Delta1]tp

[Delta1]z w(,z)

,

v w K

Kz

1

C

C

A

@ @i+1(,z)

@i(+[Delta1],z)

@v(,z) (40)

=

1

2 [notdef] [parenleftbigg]

(33)

[Delta1]tp

[Delta1]

2[notdef]

v( + [Delta1],z)

r2 [notdef]cos( + [Delta1]) cos() ,

@ @i+1([Delta1],z)@i(,z)

@v(,z) =

and assume constant velocities and mixing coefcients during the macro time step. To infer circulation patterns and mixing coefcients from the measurements of air densities and mixing ratios, the Jacobian matrix F,

F = f 1,...,f N

[parenrightbig]

= (fj,n) =

1

2r [notdef]

[Delta1]tp

[Delta1] [notdef]

@xI,j @qn

cos() cos( [Delta1])

(41)

(34)

[notdef]

1 +

[Delta1]tp

[Delta1] [notdef]

v( [Delta1],z)

r +

[Delta1]tp

[Delta1]z w( [Delta1],z)

,

,

is needed, where N = 4J , where J =

Pgases0 Jg, because there are four unknown quantities, vj , wj , K;j , and Kz;j ,

at each grid point of the nominal region where these variables shall be inferred. The elements of F are calculated from

=

@xI

@q1 ,...,

@xI @qN

@ @i+1(,z+[Delta1]z)@i(,z)

@v(,z) (42)

=

1

2 [notdef]

w(,z) r ,

www.atmos-chem-phys.net/16/14563/2016/ Atmos. Chem. Phys., 16, 1456314584, 2016

[Delta1]tp

[Delta1]z [notdef]

[Delta1]tp

[Delta1] [notdef]

r2(r + [Delta1]z)2

14572 T. von Clarmann and U. Grabowski: Circulation and mixing from tracer measurements

@ @i+1(,z)

@i(,z+[Delta1]z)

@v(,z) (43)

=

1

2 [notdef]

@ @i+1(,z)

@i(+[Delta1],z)

@w(,z) (53)

=

1

2 [notdef]

[Delta1]tp

[Delta1]z [notdef]

[Delta1]tp

[Delta1] [notdef]

(r + [Delta1]z)2

r3 w(,z + [Delta1]z),

@ @i+1(+[Delta1],z)@i(,z+[Delta1]z)@v(,z) =

[Delta1]tp

[Delta1] [notdef]

[Delta1]tp

[Delta1]z [notdef]

v( + [Delta1],z)

r [notdef]

cos( + [Delta1])

cos() ,

1

2 [notdef]

[Delta1]tp

[Delta1] [notdef]

[Delta1]tp

[Delta1]z (44)

[notdef]

@ @i+1(+[Delta1],z)

@i(+[Delta1],z[Delta1]z)

@w(,z) (54)

=

1

2 [notdef]

(r + [Delta1]z)2

r3 [notdef]

cos() cos( + [Delta1])

w(,z + [Delta1]z),

[Delta1]tp

[Delta1] [notdef]

[Delta1]tp

[Delta1]z [notdef]

v(,z)

r [notdef]

cos() cos( + [Delta1])

,

@ @i+1([Delta1],z+[Delta1]z)@i(,z)

@v(,z) =

1

2 [notdef]

[Delta1]tp

[Delta1] [notdef]

[Delta1]tp

[Delta1]z (45)

[notdef]

@ @i+1(,z+[Delta1]z)@i(,z)

@w(,z) =

1

2 [notdef]

[Delta1]tp

[Delta1]z [notdef]

r2(r + [Delta1]z)2

(55)

r

(r + [Delta1]z)2 [notdef]

cos() cos( [Delta1]) [notdef]

w( [Delta1],z),

[notdef]

1 +

[Delta1]tp

[Delta1] [notdef]

v(,z)

r + 2[Delta1]tp

[Delta1]z [notdef] w(,z)

,

@ @vmri+1(,z)@vmri(,z)

@v(,z) (46)

=

[Delta1]tp [Delta1]

@ @i+1(,z)

@i(,z+[Delta1]z)

@w(,z) (56)

=

1

2 [notdef] [parenleftbigg]

2[notdef]1 r2 [notdef]

v(,z) +

v( [Delta1],z) 2

,

[Delta1]tp

[Delta1]z

2[notdef]

(r + [Delta1]z)2r2 [notdef] w(,z + [Delta1]z),

@ @i+1(,z[Delta1]z)@i(,z)

@w(,z) =

@ @vmri+1(+[Delta1],z)@vmri(+[Delta1],z)@v(,z) =

[Delta1]tp [Delta1]

2[notdef]

v( + [Delta1],z)2r2 , (47)

1

2 [notdef]

[Delta1]tp

[Delta1]z [notdef]

r2(r [Delta1]z)2

(57)

[notdef]

1 +

[Delta1]tp

[Delta1] [notdef]

v(,z [Delta1]z)

r [Delta1]z +

[Delta1]tp

[Delta1]z [notdef] w(,z [Delta1]z)

,

@ @vmri+1(,z)

@vmri([Delta1],z)@v(,z) (48)

=

1

2r [notdef]

[Delta1]tp

[Delta1] [notdef]

1 +

[Delta1]tp

[Delta1] [notdef]

v( [Delta1],z) r

@ @i+1(+[Delta1],z[Delta1]z)@i(,z)

@w(,z) =

1

2 [notdef]

[Delta1]tp

[Delta1] [notdef]

[Delta1]tp

[Delta1]z (58)

[notdef]

,

r2(r [Delta1]z)2 [notdef]

cos() cos( + [Delta1])

v(,z [Delta1]z)

r [Delta1]z

,

@ @vmri+1(+[Delta1],z)@vmri(,z)

@v(,z) =

1

2r

[Delta1]tp [Delta1]

2 v( + [Delta1],z)

r , (49)

@ @i+1(,z+[Delta1]z)@i(+[Delta1],z)

@w() =

1

2 [notdef]

[Delta1]tp

[Delta1] [notdef]

[Delta1]tp

[Delta1]z (59)

[notdef]

@ @vmri+1(,z)

@vmri(+[Delta1],z)@v(,z) =

[Delta1]tp

[Delta1] [notdef]

1 r

12

[Delta1]tp

[Delta1] [notdef] v(,z) r

r2(r + [Delta1]z)2 [notdef]

cos( + [Delta1])

cos() [notdef]

v( + [Delta1],z)

r ,

. (50)

Entries not mentioned here are zero. Entries relevant to w are

@ @i+1(,z)@i(,z)

@w(,z) =

[Delta1]tp

[Delta1] [notdef]

[Delta1]tp

[Delta1]z [notdef]

v(,z)

r (51)

@ @vmri+1(,z)@vmri(,z)

@w(,z) (60)

=

[Delta1]tp [Delta1]z

2[notdef]

w(,z) +

w(,z [Delta1]z) 2

,

[Delta1]tp [Delta1]z

2w(,z) 1 2

[Delta1]tp [Delta1]z

2w(,z [Delta1]z),

@ @i+1(,z+[Delta1]z)@i(,z+[Delta1]z)@w(,z) =

1

2

@ @vmri+1(,z+[Delta1]z)@vmri(,z+[Delta1]z)

@w(,z) =

[Delta1]tp [Delta1]z

2[notdef]

w(,z + [Delta1]z)2 , (61)

[Delta1]tp [Delta1]z

2w(,z + [Delta1]z), (52)

Atmos. Chem. Phys., 16, 1456314584, 2016 www.atmos-chem-phys.net/16/14563/2016/

T. von Clarmann and U. Grabowski: Circulation and mixing from tracer measurements 14573

@ @vmri+1(,z)

@vmri(,z[Delta1]z)

@w(,z) (62)

=

1

2 [notdef]

@ @vmri+1(,z+[Delta1]z)@vmri(,z)

@Kz(,z) =

1

2 [notdef]

[Delta1]tp ([Delta1]z)2 [notdef]

r + [Delta1]z2

2(r + [Delta1]z)2, (73)

@ @vmri+1(,z)

@vmri(,z[Delta1]z)

@Kz(,z) =

[Delta1]tp

[Delta1]z [notdef]

1 + w(,z [Delta1]z) [notdef] [Delta1]tp [Delta1]z

,

1

2 [notdef]

[Delta1]tp ([Delta1]z)2 [notdef]

r [Delta1]z2

2r2 , (74)

@ @vmri+1(,z)

@vmri(,z+[Delta1]z)

@Kz(,z) =

@ @vmri+1(,z+[Delta1]z)@vmri(,z)

@w(,z) =

1

2 [notdef] [parenleftbigg]

[Delta1]tp

[Delta1]z

2[notdef] w(,z + [Delta1]z), (63)

1

2 [notdef]

[Delta1]tp ([Delta1]z)2 [notdef]

r + [Delta1]z2

2r2 , (75)

@ @vmri+1(,z[Delta1]z)@vmri(,z)

@Kz(,z) =

@ @vmri+1(,z)

@vmri(,z+[Delta1]z)

@w(,z) (64)

=

1

2 [notdef]

2(r [Delta1]z)2. (76)

With these derivatives we linearize the prediction with respect to wind and mixing coefcients; i.e., we linearly predict the new atmospheric state for a given initial state q0 as a function of wind and mixing ratios.

xI = x0 + F(q q0) (77)

This formulation gives access to the winds and mixing ratios via inversion of F.

Assuming linearity and Gaussian statistics, the most likely set q of winds and mixing ratios during the macro time interval minimizes the following cost function:

[notdef]21 = xm;I xI

TS1r

1

2 [notdef]

[Delta1]tp ([Delta1]z)2 [notdef]

r [Delta1]z2

[Delta1]tp

[Delta1]z [notdef]

1 2([Delta1]tp)

([Delta1]z) [notdef] w(,z)

.

Entries relevant to K are

@ @vmri+1(,z)@vmri(,z)

@K(,z) =

[Delta1]tp ([Delta1])2 [notdef]

1

2r2 (65)

[notdef]

cos

+ [Delta1]2 [parenrightBig] + cos

[Delta1]2 [parenrightBig]

cos() ,

@ @vmri+1( [Delta1],z)@vmri( [Delta1],z)@K(,z) =

[Delta1]tp ([Delta1])2 [notdef]

1

2r2 [notdef]

cos

[Delta1]2 [parenrightBig]

cos( [Delta1])

, (66)

xm;I xI

[parenrightbig]

(78)

xm;I x0 F(q q0)

TS1r

@ @vmri+1(,z)

@vmri([Delta1],z)

@K(,z) =

1

2r2 [notdef]

[Delta1]tp ([Delta1])2 [notdef]

cos( [Delta1]2 )cos() , (67)

@ @vmri+1(,z)

@vmri(+[Delta1],z)

@K(,z) =

xm;I x0 F(q q0)

[parenrightbig]

,

where xm;I is the measured state at the end of the macro

time step and Sr is the error covariance matrix of the residual, which is, under the assumption that prediction error and measurement errors are uncorrelated, the sum of the prediction covariance matrix and the measurement covariance matrix, both after the macro time step:

Sr = Sm,I + Sp, (79)

where Sp is a J [notdef] J matrix containing those elements of SI

which are relevant to xI . Sm,I is the measurement error covariance matrix of the atmospheric state after the macro time step. The minimization of the cost function gives the following estimate q of winds and mixing coefcients:

q = q0 +

FTS1rF

1

2r2 [notdef]

[Delta1]tp ([Delta1])2 [notdef]

cos( + [Delta1]2 )cos() , (68)

@ @vmri+1(+[Delta1],z)@vmri(,z)

@K(,z) =

1

2r2 [notdef]

[Delta1]tp ([Delta1])2 [notdef]

cos( + [Delta1]2 ) cos( + [Delta1])

, (69)

@ @vmri+1([Delta1],z)@vmri(,z)

@K(,z) =

1

2r2 [notdef]

[Delta1]tp ([Delta1])2 [notdef]

cos( [Delta1]2 ) cos( [Delta1])

. (70)

And nally, entries relevant to Kz are

@ @vmri+1(,z)@vmri(,z)

@Kz(,z) =

[Delta1]tp ([Delta1]z)2 [notdef]

1FTS1r(xm;I xI ). (80)

The matrix FTS1rF can be singular either because the related system of equations is under-determined or ill-posed due to nearly linearly dependent equations. Singularity is addressed by adding the following constraint term to the cost function of Eq. (78):

www.atmos-chem-phys.net/16/14563/2016/ Atmos. Chem. Phys., 16, 1456314584, 2016

r + [Delta1]z2

2+

22r2 , (71)

@ @vmri+1(,z [Delta1]z)@vmri(,z [Delta1]z)

@Kz(,z) =

r [Delta1]z2

[Delta1]tp ([Delta1]z)2 [notdef]

r [Delta1]z2

22(r [Delta1]z)2, (72)

14574 T. von Clarmann and U. Grabowski: Circulation and mixing from tracer measurements

[notdef]2con = q qa

xm;I xI,it + Fit(qit qa)

[parenrightbig]

.

(81) [notdef]2 = [notdef]21 + [notdef]2con, (82)

where qa is some prior assumption on velocities and mixing coefcients. R is a J [notdef] J regularization matrix of which the

choice is discussed below. From this, the constrained estimate of velocities and mixing coefcients can be inferred:

q = qa + FTS1rF + R[parenrightBig]1FTS1r(xm;I xI ). (83)

An equivalent formulation, which is more efcient if the dimension of q is larger than that of x (under-determined problem) but which requires a non-singular regularization matrix and does not give easy access to diagnostics (see below), is (Rodgers, 2000)

q = qa + R1FT[parenleftBig]

FR1FT + Sr,I [parenrightBig]1(xm;I xI ). (84)

The covariance matrix characterizing the uncertainty of estimated winds and mixing coefcients is

Sq = FTS1rF + R[parenrightBig]1 FTS1rF FTS1rF + R[parenrightBig]1, (85)

and the estimated winds and mixing coefcients are related to the true ones as

A =

@ q

@q =

TR

q qa [parenrightbig]

FTS1rF + R[parenrightBig]1FTS1rF, (86)

which is unity in the case of the unconstrained estimation of q. In the case of Newtonian iteration, Eqs. (85)(86) are evaluated using the Jacobian F valid at the solution.

Due to the concentration dependence of the source function and the q dependence of F, Eq. (29) is valid only in linear approximation. This is helped by putting the inversion in the context of a Newtonian iteration (see, e.g., Rodgers, 2000, p. 85). Equation (80) becomes

qit+1 = qit +

FTitS1rFit

With qa = 0 and diagonal R = I we get the smallest pos

sible velocities and mixing coefcients still consistent with the measurement, where tuning parameter will be set depending on how large t residuals may be for the user to still consider them to be consistent. With R being diagonally blockwise composed of squared and scaled rst-order nite difference operators and qa = 0, smooth elds of wind

vectors and mixing coefcients can be enforced. Setting qa the result of the previous macro time step corresponds to sequential data assimilation. In this application R is set to the reciprocal uncertainty of qa plus some margin to allow for variability of velocity and mixing coefcients in time. And nally, if prior knowledge is formed by independent measurements and their reciprocal uncertainties as a constraint matrix, or within the debatable framework of Bayesian statistics, estimates q would even be the most probable estimate of

velocities and mixing ratios.

5 Proof of concept

5.1 Prediction of the atmospheric state

In a rst step we test the predictive power of the formalism dened by Eqs. (3)(29). Since the formalism itself is deductive and starts from a well-established theoretical concept, the purpose of the test is solely to verify that the implementation of the formalism is correct and that involved numerical approximations are adequate. As a consequence of the Bonini paradox (see Bonini, 1963, and Starbuck, 1975), a model is the harder to understand the more complex it is.While the predictive power of a model usually increases with complexity, this does not necessarily hold for its explanatory power. Thus, we have decided to test our model on the basis of very simple test cases, where major failure of the model is immediately obvious. Four test cases have been chosen, each dedicated to one kinematic variable (v, w, K and Kz), while the other three were set to zero.

In the rst case, v was set to 1/10 of the Courant limit (Courant et al., 1928) (about 0.17/cos() m s1) everywhere. As one would expect from the continuity equation applied to a spherical atmosphere, no changes in air density except boundary effects at the poles were observed, and structures near the equator were transported by about 4 within a month, as expected from the equation of motion. A Gaussian-shaped perturbation of a half width of one latitudinal grid width (4 ) causes an upwind wiggle of less than 0.7 % of the amplitude of the perturbation at a meridional velocity of one grid point per month. There is no discernable change in the width of the transported structure. Similarly, for the second case a constant eld of w of 1.1 [notdef] 103 m s1 lifts a struc

ture upwards by about 3 km per month. Mixing coefcients were veried to smear out structures in the respective direc-

Atmos. Chem. Phys., 16, 1456314584, 2016 www.atmos-chem-phys.net/16/14563/2016/

1FTitS1r(xm;I xI,it), (87)

where it is the iteration index. Equation (83) becomes

qit+1 = qit + FTitS1rFit + R[parenrightBig]1

(88)

FTitS1r(xm;I xI,it) R(qit qa)[parenrightBig] or alternatively

qit+1 = qa +

FTitS1rFit + R[parenrightBig]1

(89)

FTitS1r xm;I xI,it + Fit(qit qa)

[parenrightbig]

,

and Eq. (84) becomes

qit+1 = qa + R1FTit[parenleftBig]

FitR1FTit + Sr,I [parenrightBig]1

(90)

T. von Clarmann and U. Grabowski: Circulation and mixing from tracer measurements 14575

tion while leaving air density and structures in the orthogonal direction unchanged.

5.2 Inversion of simulated measurements

Case studies based on real measurement data are inadequate as the sole proof of concept because the truth is unknown and the result thus is unveriable. Instead we rst test our scheme on the basis of simulated atmospheric states and consider the scheme as veried if the velocities and mixing coefcients used to simulate the atmospheric states are sufciently well reproduced. In the noise-free well-conditioned case, one might even expect, within the numerical precision of the system, the exact reproduction of the reference data; due to the weak but non-zero dispersivity of the numerical transport scheme, the wiggles discussed in the previous subsection cause, at some grid points, D matrix entries of the wrong sign. In order to ght resulting convergence problems of the inversion, some small regularization is at least necessary, even if the system of equations to be solved is well or overdetermined. Since the system in reality tends to be ill-conditioned and the constraint applied to the inversion prevents the reproduction of the reference data, we use a variety of idealized tracers instead. After this initial test of functionality, more and more realistic test cases are constructed in order to study the competing inuence of constraint and measurement data on the solution.

During code development, a series of basic tests of increasing complexity were performed, including a variety of mixing ratio distributions transported with various velocity elds. The main lesson learned was that, even when the rows of Eq. (80) were independent and no formal ill-posedness could be diagnosed, the unphysical upwind wiggles in the vicinity of the mixing ratio peak as discussed in the previous subsection could trigger errors which are boosted during the iteration. This problem, which is associated with sharp structures and large velocities (of the order of one grid width per macro time step) can be solved by the use of a smoothing regularization matrix R as discussed in the last paragraph of

Sect. 4, however at the cost of degraded spatial resolution of the result.

As an example we show the following test. An altitude-independent meridional velocity eld (in degrees per month)

v =

23.6562

2 cos (91)

was chosen, while the vertical velocity was set to zero (Fig. 2a). This particular choice of the meridional velocity eld keeps boundary problems at the poles due to divergence in a circulation which is not closed reasonably small. Initial mixing ratio distributions (in ppmv) of four idealized gases were constructed such that gas (a) was sugarloaf shaped, while the mixing ratio distributions of the gases (b) to (d) were monotonic slopes.

vmra = e[parenleftBig]

z36

1.5

2/72

+5.0 [notdef] 104[parenrightBigg] e(

6(+6)

40 )2/12, (92)

(30 + 0.02z2 0.2)3

1000 , (93) vmrc = 20 [notdef] e0.02z+0.03, (94)

and

vmrd = 30 [notdef] e0.03z0.002, (95)

where z is altitude above surface in km and is latitude in degrees. The distribution of gas (a) is shown in Fig. 2b. The transport problem was solved in the forward mode using the velocity eld as dened in Eq. (91) for integrating the tendency equation (Eq. 28) over 1 month. The resulting distribution of gas (a) after this macro time step is shown Fig. 2c.The sugarloaf is transported to the north by the expected distance. Due to the latitudinal gradients with lower velocities in the luff of the sugarloaf and larger velocities in the lee, the shape of the sugarloaf is slightly distorted, leading to a steeper slope in the luff and a atter slope in the lee. No indication of problems with diffusivity or dispersivity of the transport scheme is seen. At the poles boundary problems are visible which are unavoidable when such an unrealistic (but instructive) velocity eld is used. The inversion nicely reproduces the initial velocity eld (not shown because hardly discernable from Fig. 2a). Due to the smoothing constraint the peak velocity is decreased by 18 %. This smoothing is the price to pay for the stabilization of the retrieval in the presence of the boundary problems discussed above. The residual between the simulated measurement of the distribution of gas(a) at t0 +[Delta1]t and the distribution predicted with the resulting

velocity eld is shown in (Fig. 2d).

5.3 Case study with MIPAS measurements

The risk of case studies based on simulated data typically is that not all difculties encountered with real data are foreseen during theoretical studies. In order to demonstrate applicability to real data, global monthly latitudealtitude distributions of CFC-12, CH4, N2O, and SF6 (Kellmann et al., 2012; Plieninger et al., 2015; Haenel et al., 2015) measured with the MIPAS (Fischer et al., 2008) were used. The purpose of these tests is to demonstrate the feasibility of the method presented. An investigation of the atmospheric circulation on the basis of this method applied to MIPAS data is left for a companion paper. For this proof of concept, sinks of these long-lived tracers have been ignored, but these will certainly be considered in scientic applications.

For this case study, zonal monthly mean distributions of air densities and mixing ratios of these four species from September and October 2010 were used. Figure 3 shows the measured distributions of these quantities in September (left

www.atmos-chem-phys.net/16/14563/2016/ Atmos. Chem. Phys., 16, 1456314584, 2016

vmrb =

14576 T. von Clarmann and U. Grabowski: Circulation and mixing from tracer measurements

Vj,Vz

0 -90 -60 -30 0 30 60 90Latitude/deg

column) and October (middle column) and the residuals between the measured and predicted contributions for October (right column). The residuals are reasonably small and show, except for methane in the polar upper stratosphere, no patterns which would suggest peculiarities with the inferred kinematic quantities.

The resulting circulation vectors which best explain the change in the mixing ratio distributions from September to October 2010 are shown in the upper left panel of Fig. 4.Winter polar subsidence, summer polar upwelling, the mesospheric overturning circulation, the upper and lower branches of the BrewerDobson circulation, and the tropical pipe are clearly visible. Details of the tropical pipe are visible in the right upper panel. As expected, the BrewerDobson circulation is much more pronounced in the northern (early) winter hemisphere. Velocities are roughly consistent with the mean ages of stratospheric air as determined by Stiller et al. (2008) and Haenel et al. (2015) in the sense that the quotient of the typical circulation velocity and the distance between equator and pole gives an age estimate of the correct order of magnitude. While the inferred eld of circulation vectors shows many detailed features demanding scientic investigation in their own right (e.g., the latitude offset between the intertropical convergence zone and the stratospheric tropical pipe or the interface between the stratospheric two-cell circulation and the overturning mesospheric circulation and the transition altitude between them), the reproduction of the expected features justies condence in the method proposed. Resulting mixing coefcients K and Kz are shown in the left and right lower panels, respectively. Negative mixing coefcients indicate counter-gradient mixing, which seems to be most pronounced in the tropical upper stratosphere.

Jacobian elements with respect to v values and K values seem to form a null space. Thus the K values were constrained to zero using diagonal components only in the R matrix diagonal blocks associated with the K values. The strength of this constraint was adjusted such that the K values were as small as possible as long as this did not boost the residual. Resulting K and Kz distributions are shown in Fig. 4.

The errors in the estimated transport velocities and mixing coefcients have been estimated according to Eq. (85) and are shown in Fig. 5. The uncertainties in the transport velocities caused by the propagation of measurement errors are in the 1 % range, indicating that the information contained in the measurements is adequate for the purpose of retrieving circulation parameters. It even seems possible to improve the time resolution of the circulation analysis and aim at weekly instead of monthly temporal sampling.

Larger errors above 65 km altitude and at the bins closest to the pole are border effects, resulting from the fact that no symmetric derivatives can be calculated there. The uncertainties in K show the same patterns as the K values themselves.

Atmos. Chem. Phys., 16, 1456314584, 2016 www.atmos-chem-phys.net/16/14563/2016/

80

max|Vj|=11.82 degmth , max|Vz|= 0.00 kmmth

70

(a)

10

60

Altitude/km

50

8

40

6

30

4

20

10

2

vmr_a/ppmv, t=t

-90 -60 -30 0 30 60 90 Latitude/deg

60

(b)

50

0.5

Altitude/km

40

0.0

30

20

-0.5

10

-1.0

vmr_a/ppmv, t=t + Dt

-90 -60 -30 0 30 60 90 Latitude/deg

60

(c)

50

0.5

Altitude/km

40

0.0

30

20

-0.5

10

-1.0

Residuum vmr_a/ppmv

-90 -60 -30 0 30 60 90 Latitude/deg

60

(d)

0.2

0.0

50

-0.2

Altitude/km

40

-0.4

30

-0.6

20

-0.8

-1.0

10

-1.2

Figure 2. Test case: (a) velocity eld; (b) initial distribution of gas (a); (c) distribution of gas (a) after 1 month; (d) residual between distribution of gas (a) after 1 month predicted with the resulting and the correct velocities. The units of the color bar in (a) are degrees per month.

T. von Clarmann and U. Grabowski: Circulation and mixing from tracer measurements 14577

Air density/ISA76

-90 -60 -30 0 30 60 90 Latitude/deg

Air density/ISA76

-90 -60 -30 0 30 60 90 Latitude/deg

Residuum air density kg m

-90 -60 -30 0 30 60 90 Latitude/deg

60

60

60

50

50

50

Altitude/km

altitude/km

altitude/km

40

40

40

30

30

30

20

20

20

10

10

10

CFC12/ppbv

-90 -60 -30 0 30 60 90 Latitude/deg

CFC12/ppbv

-90 -60 -30 0 30 60 90 Latitude/deg

Residuum CFC12/ppbv

-90 -60 -30 0 30 60 90 Latitude/deg

60

60

60

50

50

50

Altitude/km

altitude/km

altitude/km

40

40

40

30

30

30

20

20

20

10

10

10

CH /ppmv

-90 -60 -30 0 30 60 90 Latitude/deg

CH /ppmv

-90 -60 -30 0 30 60 90 Latitude/deg

Residuum CH /ppmv

-90 -60 -30 0 30 60 90 Latitude/deg

60

60

60

50

50

50

Altitude/km

altitude/km

altitude/km

40

40

40

30

30

30

20

20

20

10

10

10

N O/ppmv

-90 -60 -30 0 30 60 90 Latitude/deg

N O/ppmv

-90 -60 -30 0 30 60 90 Latitude/deg

Residuum N O/ppbv

-90 -60 -30 0 30 60 90 Latitude/deg

60

60

60

50

50

50

Altitude/km

altitude/km

altitude/km

40

40

40

30

30

30

20

20

20

10

10

10

SF /pptv

-90 -60 -30 0 30 60 90 Latitude/deg

SF /pptv

-90 -60 -30 0 30 60 90 Latitude/deg

Residuum SF /pptv

-90 -60 -30 0 30 60 90 Latitude/deg

60

60

60

50

50

50

Altitude/km

altitude/km

altitude/km

40

40

40

30

30

30

20

20

20

10

10

10

Figure 3. Measured distributions in September (left column) and October (middle column) and residual distributions between October measurements and predictions for October (right column) for air density and mixing ratios of CFC-12, CH4, N2O, and SF6 (top to bottom).

Grey grid boxes indicate nonavailability of valid data.

www.atmos-chem-phys.net/16/14563/2016/ Atmos. Chem. Phys., 16, 1456314584, 2016

14578 T. von Clarmann and U. Grabowski: Circulation and mixing from tracer measurements

step by measured value, depletion does not accumulate, even if no sink functions are considered) and

d. the circular reasoning that the lifetimes of nonideal tracers depend on their trajectories (and thus atmospheric circulation), while the determination of the circulation requires knowledge of the lifetimes, can be resolved; our scheme requires knowledge only of the local, not the global, lifetimes

e. the method is an empirical method which does not involve any dynamical model; i.e., the forces which cause the circulation are not required.

The method only nds that kinematic state of the atmosphere which, according to the continuity equation, t the measurements best. These kinematic state values are provided as model diagnostics to assess the performance of dynamical models. Due to these advantages, the major problems in the empirical analysis of the BrewerDobson circulation as mentioned by Butchart (2014) are solved. Problems related to our method are

a. sensitivity of the inferred kinematic quantities to locally varying biases,

Atmos. Chem. Phys., 16, 1456314584, 2016 www.atmos-chem-phys.net/16/14563/2016/

Figure 4. Resulting circulation vectors (v(z,),w(z,)) (upper left panel), where colors on the red side of the rainbow color scale represent higher velocities; a zoomed-in view of this (upper right panel); mixing coefcients K (lower left panel) and Kz (lower right panel).

6 Discussion

The analysis of the age of stratospheric air can be understood as an integrated look at the equations of motion of stratospheric air because the total travel time of the air parcel through the stratosphere is represented. The renement of this method which analyzes the mean age just considers a weighted mean of the above, but it is still an integral method. Contrary to these integral methods, our direct inversion scheme supports a in approximation, due to discrete sampling in the time domain differential view of the same problem. The related advantages are the following:

a. independence of assumptions on the age spectrum because during each time step mixing is explicitly considered

b. insensitivity to SF6 depletion in the mesosphere (see, e.g., Reddmann et al., 2001; Stiller et al., 2012) because the scheme uses the actual entry values of subsiding air as a reference

c. applicability to nonideal tracers in the stratosphere (since the atmospheric state is updated for each time

T. von Clarmann and U. Grabowski: Circulation and mixing from tracer measurements 14579

SD Vj/10 degmth

-90 -60 -30 0 30 60 90 Latitude/deg

SD Vz/10 kmmth

-90 -60 -30 0 30 60 90 Latitude/deg

60

60

60

80

50

50

50

60

Altitude/km

Altitude/km

40

40

40

30

40

30

30

20

20

20

20

10

10

10

SD Kj/10 m s

-90 -60 -30 0 30 60 90 Latitude/deg

SD Kz/10 m s

-90 -60 -30 0 30 60 90 Latitude/deg

3.5

60

60

8

3.0

50

50

2.5

6

Altitude/km

Altitude/km

40

40

2.0

30

4

30

1.5

1.0

20

20

2

0.5

10

10

Figure 5. Estimated uncertainties of v(z,) (upper left panel), w(z,) (upper right panel), K (lower left panel), and Kz (lower right panel).

b. a tendency towards ill-posedness of the inversion if distributions of too few tracers with too similar a morphology are used, and

c. the usual artifacts arising if the numerical discretization chosen is too coarse.

Results of the case study presented in Sect. 5.3 suggest that these problems are under control in the current application of the proposed scheme.

7 Conclusions and outlook

We have presented a method which infers mixing coefcients and effective velocities of a 2-D atmosphere by inversion of the continuity equation. The main steps of this procedure are

a. the integration of the continuity equation over time to predict pressure and mixing ratios for given initial pressures and mixing ratios and initially guessed velocities and mixing coefcients;

b. the propagation of errors of initial pressures and mixing ratios onto the predicted pressures and mixing ratios, by differentiation of the predicted state with respect to the initial state and generalized Gaussian error propagation;

c. the estimation of the sensitivities of the predicted state with respect to the velocities and mixing coefcients; and

d. the minimization of a quadratic cost function involving the residual between measured and predicted state at the end of the forecasting interval by inversion of the continuity equation.

The inferred velocities are suggested to be used as a model diagnostic in order to avoid problems encountered with other model diagnostics like mean age of stratospheric air. It is important to note that the diagnostics inferred here are effective transport velocities and effective mixing coefcients in the sense that they include eddy transport and diffusion terms.Thus, they cannot simply be compared to zonal mean velocities and mixing coefcients of a 3-D model, but the eddy terms have to be considered when these diagnostics quan-

www.atmos-chem-phys.net/16/14563/2016/ Atmos. Chem. Phys., 16, 1456314584, 2016

14580 T. von Clarmann and U. Grabowski: Circulation and mixing from tracer measurements

tities are calculated. The application of this method to SF6 distributions measured by MIPAS (Stiller et al., 2012) to diagnose the BrewerDobson circulation is discussed in a companion paper. Obvious future activities are the extension of this method to three dimensions and the inclusion of sink functions of non-inert species to explore a larger number of tracers in order to better constrain the related inverse problem.

8 Data/Code availability

The code is still under development. The data used are compiled in Supplement 2. The complete MIPAS data are available at http://www.imk-asf.kit.edu/english/308.php

Web End =http://www.imk-asf.kit.edu/english/308.php http://www.imk-asf.kit.edu/english/308.php

Web End = .

Atmos. Chem. Phys., 16, 1456314584, 2016 www.atmos-chem-phys.net/16/14563/2016/

T. von Clarmann and U. Grabowski: Circulation and mixing from tracer measurements 14581

Appendix A: From 3-D to 2-D

The reduction of the transport problem from three to two dimensions involves Reynolds decomposition of the three-dimensional continuity equation and subsequent zonal averaging and gives rise to eddy mixing and transport terms.The inference of effective two-dimensional transport velocities and effective mixing coefcients from measurements discussed in the main paper relies on the fact that, within certain assumptions and approximations, all these eddy effects can be understood as additional pseudo-advection and pseudo-mixing terms according to the advection equation and Ficks law, with gas-independent pseudo-velocities and pseudo-mixing coefcients. The exact interpretation of the two-dimensional velocities and mixing coefcients inferred from the measurements depends on the approximations made.We apply our scheme to zonally averaged mixing ratios (no mass-weighted averaging). Contrary to the main text, where the symbols v, w, K, and Kz are used in the two-dimensional system, here zonal averages are indicated by a bar.

Assuming

that the deviations from the zonal mean are small compared to the zonal mean itself such that linearization is justiable,

that meridional advection is negligibly small compared to zonal advection,

that the time variation of the zonal mean quantities is assumed to be much slower than the time variation of the deviations from the zonal mean, which corresponds to the assumption of a quasi-steady state,

Tung (1982) derives a two-dimensional approximation to the continuity equation which, adjusted to our notation, written for geometric altitudes instead of potential temperatures, and using the shallowness approximation, reads

@

@t vmrg +

v

@

@t [prime] [prime]

[parenrightbigg]1r

@

@ vmrg (A1)

+

w

@

@t [prime][Phi1][prime] [parenrightbigg]

@ @zvmrg

[parenleftbigg]

@

@t +

u

r cos

@ @

[Phi1][prime] = w[prime], (A3)

and

[parenleftbigg]

@

@t +

u

r cos

@ @

[prime] = S[prime]. (A4)

Further,

K =

1 (v)[prime] [prime], (A5)

Kzz =

1 (w)[prime][Phi1][prime], (A6)

Kz =

1 (v)[prime][Phi1][prime], (A7)

and

Kz =

1 (w)[prime] [prime]. (A8)

Assuming further that

our scheme is applied only to long-lived species, such that chemical eddy terms can be ignored because chemical lifetimes are long compared to transport lifetimes (see Pyle and Rogers, 1980) and

wave disturbances are dominated by steady or periodic terms, such that the terms with the mixed second derivative terms tend to disappear (Matsuno, 1980; Clark and Rogers, 1978; Pyle and Rogers, 1980; Tung, 1982) and an even weaker approximation is sufcient, namely

@

@t [prime][Phi1][prime] 0, (A9)

making it equivalent with the assumption that Kz

Kz,

Eq. (A1) can be rewritten as a tendency equation of the type

@

@t vmrg =

v r

@

@ vmrg w

@

@zvmrg (A10)

+

1 r2

@ @

K

@

@ vmrg [bracketrightbigg]

+

1 r

@ @

K 1r

@

@ vmrg

[parenrightbigg]

@ @z

K z

@ @zvmrg

,

1 r

@ @

Kz @ @zvmrg

@ @z

Kz 1r

@

@ vmrg [parenrightbigg]

where

v = v

1

@

@t [prime] [prime]

1 r

@

@ K (A11)

@ @z

Kzz @ @zvmrg

[parenrightbigg]

@

@z Kz,

@ @zKz

1

= S

1 r

@

@ ((v)[prime][prime])

@ @z((z)[prime][prime])

@

@t ([prime][prime]),

@

@z Kzz (A12)

w = w

1

@

@t [prime][Phi1][prime]

1

where [prime], [Phi1][prime] and [prime] are dened by

[parenleftbigg]

@

@t +

@

@ Kz,

www.atmos-chem-phys.net/16/14563/2016/ Atmos. Chem. Phys., 16, 1456314584, 2016

u

r cos

@ @

1 r

@

@ Kz

1 r

[prime] = v[prime], (A2)

14582 T. von Clarmann and U. Grabowski: Circulation and mixing from tracer measurements

K = K, (A13)

and

K z = Kzz (A14)

can be understood as virtual velocities and virtual mixing coefcients. All terms in Eq. (A10) have the mathematical structure of either an advection equation or a Ficks law.Thus, Eqs. (A11)(A14) provide the interpretation of the velocities and mixing coefcients inferred in the main part of the paper. These can be identied with the virtual velocities and virtual mixing coefcients inferred above. Obviously, this interpretation depends on the approximations made, and different approximations would lead to a different interpretation. If mixing coefcients describing subscale effects in the original 3-D model (the last two terms in Eq. 3) are to be considered, then the effective mixing coefcient is the sum of the mixing coefcients describing these subscale effects plus the respective K term accounting for the eddy mixing.

We would like to emphasize that none of the approximations and assumptions discussed above are used in our proposed method to infer velocities from zonal mean mixing

ratio measurements. The discussion in this Appendix only tries to relate the resulting velocities to the velocities in a 3-D world. The ambiguities in the interpretation of the inferred effective 2-D velocities suggests that it might be promising to switch from a theoretical to an empiricist view and to no longer conceive of the zonal mean of the 3-D velocities as the true 2-D velocities but as those which satisfy the 2-D continuity equation. These can be admittedly indirectly observed with our suggested method and can be used to validate 2-D models, including their underlying concept of solving the 2-D transport problem. With this, the adequacy of the assumptions made to approximate away the headache terms which can be expressed neither as advection nor as Ficks law terms can also be tested by means of a comparison of the measured and 2-D modeled effective, i.e., transport-relevant, velocities. This empiricist turn in argumentation might not fully solve all aspects of the problem of interpretation of the observed 2-D velocities from a 3-D perspective, but at least it moves the problem from the desk of the observation scientist onto the desk of the 2-D modeler.

Atmos. Chem. Phys., 16, 1456314584, 2016 www.atmos-chem-phys.net/16/14563/2016/

T. von Clarmann and U. Grabowski: Circulation and mixing from tracer measurements 14583

The Supplement related to this article is available online at http://dx.doi.org/10.5194/acp-16-14563-2016-supplement

Web End =doi:10.5194/acp-16-14563-2016-supplement .

Acknowledgements. The authors thank two anonymous reviewers for their thorough examination of the original manuscript and for helpful and important comments. Furthermore, T. von Clarmann wishes to thank Hendrik Elbern, Richard Menard, Peter Braesicke, Bjrn-Martin Sinnhuber, and Thomas Birner for drawing his attention to some important literature and for encouragement as well as Arne Babenhauserheide for helpful discussions.

The article processing charges for this open-access publication were covered by a ResearchCentre of the Helmholtz Association.

Edited by: J.-U. Groo

Reviewed by: two anonymous referees

References

Bonini, C. P.: Simulation of information and decision systems in the rm, Prentice-Hall, Englewood Cliffs, N.J., 1963.

Brasseur, G. and Solomon, S.: Aeronomy of the Middle AtmosphereChemistry and Physics of the Stratosphere and Mesosphere, Atmospheric and Oceanographic Sciences Library 32, 3rd Edn., Springer, the Netherlands, 2005.

Butchart, N.: The Brewer-Dobson Circulation, Rev. Geophys., 52,157184, doi:http://dx.doi.org/10.1002/2013RG000448

Web End =10.1002/2013RG000448 http://dx.doi.org/10.1002/2013RG000448

Web End = , 2014.

Butchart, N., Scaife, A. A., Bourqui, M., de Grandpre, J., Hare, S.H. E., Kettleborough, J., Langematz, U., Manzini, E., Sassi, F., Shibata, K., Shindell, D., and Sigmond, M.: Simulations of anthropogenic change in the strength of the Brewer-Dobson circulation, Clim. Dynam., 27, 727741, doi:http://dx.doi.org/10.1007/s00382-006-0162-4

Web End =10.1007/s00382-006- http://dx.doi.org/10.1007/s00382-006-0162-4

Web End =0162-4 , 2006.

Chang, C. A. and St.-Maurice, J.-P.: Two-dimensional high-latitude thermospheric modeling: A comparison between moderate and extremely disturbed conditions, Can. J. Phys., 69, 10071031, doi:http://dx.doi.org/10.1139/p91-159

Web End =10.1139/p91-159 http://dx.doi.org/10.1139/p91-159

Web End = , 1991.

Clark, J. H. E. and Rogers, T. G.: The transport of trace gases by planetary waves, J. Atmos. Sci., 35, 22322235, 1978. Courant, R., Friedrichs, K., and Lewy, H.: ber die partiellen Differentialgleichungen der mathematischen Physik, Math. Ann., 100, 3274, doi:http://dx.doi.org/10.1007/BF01448839

Web End =10.1007/BF01448839 http://dx.doi.org/10.1007/BF01448839

Web End = , 1928.

Courant, R., Isaacson, E., and Rees, M.: On the Solution of Nonlinear Hyperbolic Differential Equations by Finite Differences, Comm. Pure Appl. Math., 5, 243255, doi:http://dx.doi.org/10.1002/cpa.3160050303

Web End =10.1002/cpa.3160050303 http://dx.doi.org/10.1002/cpa.3160050303

Web End = , 1952.

Engel, A., Mbius, T., Bnisch, H., Schmidt, U., Heinz, R., Levin,I., Atlas, E., Aoki, S., Nakazawa, T., Sugawara, S., Moore, F., Hurst, D., Elkins, J., Schaufer, S., Andrews, A., and Boering,K.: Age of stratospheric air unchanged within uncertainties over the past 30 years, Nat. Geosci., 2, 2831, doi:http://dx.doi.org/10.1038/ngeo388

Web End =10.1038/ngeo388 http://dx.doi.org/10.1038/ngeo388

Web End = , 2009.

Fischer, H., Birk, M., Blom, C., Carli, B., Carlotti, M., von Clar-mann, T., Delbouille, L., Dudhia, A., Ehhalt, D., Endemann, M.,

Flaud, J. M., Gessner, R., Kleinert, A., Koopman, R., Langen, J., Lpez-Puertas, M., Mosner, P., Nett, H., Oelhaf, H., Perron, G., Remedios, J., Ridol, M., Stiller, G., and Zander, R.: MIPAS: an instrument for atmospheric and climate research, Atmos. Chem. Phys., 8, 21512188, doi:http://dx.doi.org/10.5194/acp-8-2151-2008

Web End =10.5194/acp-8-2151-2008 http://dx.doi.org/10.5194/acp-8-2151-2008

Web End = , 2008. Fueglistaler, S., Dessler, A. E., Dunkerton, T. J., Folkins, I., Fu, Q., and Mote, P. W.: Tropical Tropopause Layer, Rev. Geophys., 47, RG1004, doi:http://dx.doi.org/10.1029/2008RG000267

Web End =10.1029/2008RG000267 http://dx.doi.org/10.1029/2008RG000267

Web End = , 2009.

Garcia, R. R., Randel, W. J., and Kinnison, D. E.: On the determination of age of air trends from atmospheric trace species, J. Atmos. Sci., 68, 139154, 2011.

Ghil, M.: Advances in sequential estimation for atmospheric and oceanic ows, J. Meteorol. Soc. Jpn., 75, 289304, 1997.

Ghil, M. and Malanotte-Rizzoli, P.: Data assimilation in meteorology and oceanography, Adv. Geophys., 33, 141266, 1991. Haenel, F. J., Stiller, G. P., von Clarmann, T., Funke, B., Eckert,E., Glatthor, N., Grabowski, U., Kellmann, S., Kiefer, M., Linden, A., and Reddmann, T.: Reassessment of MIPAS age of air trends and variability, Atmos. Chem. Phys., 15, 1316113176, doi:http://dx.doi.org/10.5194/acp-15-13161-2015

Web End =10.5194/acp-15-13161-2015 http://dx.doi.org/10.5194/acp-15-13161-2015

Web End = , 2015.

Hall, T. M. and Plumb, R. A.: Age as a diagnostic of stratospheric transport, J. Geophys. Res., 99, 10591070, 1994.

Hall, T. M. and Waugh, D. W.: Inuence of nonlocal chemistry on tracer distributions: Inferring the mean age of air from SF6,

J. Geophys. Res., 103, 1332713336, doi:http://dx.doi.org/10.1029/98JD00170

Web End =10.1029/98JD00170 http://dx.doi.org/10.1029/98JD00170

Web End = , 1998.

Hinkelmann, K.: Primitive Equations, in: Lectures on Numerical

Shortrange Weather Prediction, Vol. 297 of Regional Training Seminar, 306375, WMO, 1951.

Holton, J. R. and Choi, W.-K.: Transport Circulation Deduced from

SAMS Trace Species Data, J. Atmos. Sci., 45, 19291939, 1988. Jones, A. R., Thomson, D. J., Hort, M., and Devenish, B.: The U.K.

Met Ofces next-generation atmospheric dispersion model, NAME III, in: Proceedings of the 27th NATO/CCMS International Technical Meeting on Air Pollution Modelling and its Application, edited by: Borrego, C. and Norman, A.-L., Vol. XVII of Air Pollution Modeling and its Application, 580589, Springer, 2007.

Kasahara, A.: Computational Aspects of Numerical Models for

Weather Prediction and Climate Simulation, in: General Circulation Models of the Atmosphere, edited by: Chang, J., 2108, Academic Press, New York, 1977.

Kellmann, S., von Clarmann, T., Stiller, G. P., Eckert, E., Glatthor,N., Hpfner, M., Kiefer, M., Orphal, J., Funke, B., Grabowski,U., Linden, A., Dutton, G. S., and Elkins, J. W.: Global CFC-11 (CCl3F) and CFC-12 (CCl2F2) measurements with the Michelson Interferometer for Passive Atmospheric Sounding (MIPAS): retrieval, climatologies and trends, Atmos. Chem. Phys., 12, 1185711875, doi:http://dx.doi.org/10.5194/acp-12-11857-2012

Web End =10.5194/acp-12-11857-2012 http://dx.doi.org/10.5194/acp-12-11857-2012

Web End = , 2012. Krueger, A., Stremme, W., Harig, R., and Grutter, M.: Volcanic SO2 and SiF4 visualization using 2-D thermal emission spectroscopy Part 2: Wind propagation and emission rates, Atmos. Meas.

Tech., 6, 4761, doi:http://dx.doi.org/10.5194/amt-6-47-2013

Web End =10.5194/amt-6-47-2013 http://dx.doi.org/10.5194/amt-6-47-2013

Web End = , 2013. MacCormack, R. W.: The effect of viscosity in hypervelocity impact cratering, in: AIAA Hypervelocity Impact Conference, Cincinnati, Chio/30 April2 May, 1969, 69354, AIAA, American Institute of Aeronautics and Astrophysics, 1969.

www.atmos-chem-phys.net/16/14563/2016/ Atmos. Chem. Phys., 16, 1456314584, 2016

14584 T. von Clarmann and U. Grabowski: Circulation and mixing from tracer measurements

Matsuno, T.: Lagrangian motion of air parcels in the stratosphere in the presence of planetary waves, Pure Appl. Geophys., 118, 189216, 1980.

Perrin, A. and Hu, H. H.: An explicit nite-difference scheme for simulation of moving particles, J. Comput. Phys., 212, 166187, doi:http://dx.doi.org/10.1016/j.jcp.2005.06.021

Web End =10.1016/j.jcp.2005.06.021 http://dx.doi.org/10.1016/j.jcp.2005.06.021

Web End = , 2006.

Phillips, N. A.: The equations of motion for a shallow rotating atmosphere and the traditional approximation, J. Atmos. Sci., 23, 626628, 1966.

Plieninger, J., von Clarmann, T., Stiller, G. P., Grabowski, U., Glatthor, N., Kellmann, S., Linden, A., Haenel, F., Kiefer, M., Hpfner, M., Laeng, A., and Lossow, S.: Methane and nitrous oxide retrievals from MIPAS-ENVISAT, Atmos. Meas. Tech., 8, 46574670, doi:http://dx.doi.org/10.5194/amt-8-4657-2015

Web End =10.5194/amt-8-4657-2015 http://dx.doi.org/10.5194/amt-8-4657-2015

Web End = , 2015.

Pyle, J. A. and Rogers, C. F.: A modied diabatic circulation model for stratospheric tracer transport, Nature, 287, 711714, doi:http://dx.doi.org/10.1038/287711a0

Web End =10.1038/287711a0 http://dx.doi.org/10.1038/287711a0

Web End = , 1980.

Reddmann, T., Ruhnke, R., and Kouker, W.: Threedimensional model simulations of SF6 with mesospheric chemistry, J.

Geophys. Res., 106, 1452514537, doi:http://dx.doi.org/10.1029/2000JD900700

Web End =10.1029/2000JD900700 http://dx.doi.org/10.1029/2000JD900700

Web End = , 2001.

Rodgers, C. D.: Inverse Methods for Atmospheric Sounding: Theory and Practice, Vol. 2 of Series on Atmospheric, Oceanic and Planetary Physics, edited by: Taylor, F. W., World Scientic, Singapore, New Jersey, London, Hong Kong, 2000.

Salby, M. L. and Juckes, M. N.: An Algorithm for Retrieving Atmospheric Motion From Satellite Measurements of Tracer Behavior,J. Geophys. Res., 99, 14031417, 1994.

Schmidt, U. and Khedim, A.: In situ measurements of carbon dioxide in the winter Arctic vortex and at midlatitudes: An indicator of the age of stratospheric air, Geophys. Res. Lett., 18, 763 766, 1991.

Solomon, S.: Antarctic ozone: Progress towards a quantitative understanding, Nature, 347, 347354, 1990.

Starbuck, W. H.: Organizations and their environments, in: Handbook of industrial and organizational psychology, edited by: Dunnette, M. D., 10691123, Rand, Chicago, 1975.

Stiller, G. P., von Clarmann, T., Hpfner, M., Glatthor, N., Grabowski, U., Kellmann, S., Kleinert, A., Linden, A., Milz,M., Reddmann, T., Steck, T., Fischer, H., Funke, B., Lpez-Puertas, M., and Engel, A.: Global distribution of mean age of stratospheric air from MIPAS SF6 measurements, Atmos. Chem. Phys., 8, 677695, doi:http://dx.doi.org/10.5194/acp-8-677-2008

Web End =10.5194/acp-8-677-2008 http://dx.doi.org/10.5194/acp-8-677-2008

Web End = , 2008.

Stiller, G. P., von Clarmann, T., Haenel, F., Funke, B., Glatthor,N., Grabowski, U., Kellmann, S., Kiefer, M., Linden, A., Lossow, S., and Lpez-Puertas, M.: Observed temporal evolution of global mean age of stratospheric air for the 2002 to 2010 period, Atmos. Chem. Phys., 12, 33113331, doi:http://dx.doi.org/10.5194/acp-12-3311-2012

Web End =10.5194/acp-12-3311- http://dx.doi.org/10.5194/acp-12-3311-2012

Web End =2012 , 2012.

Thompson, P.: Reduction of analysis error through constraints of dynamical consistency, J. Appl. Meteorol., 8, 739742, 1969. Tung, K. K.: On the Two-Dimensional Transport of Stratspheric

Trace Gases in Isentropic Coordinates, J. Atmos. Sci., 39, 2330 2355, 1982.von Clarmann, T., Hpfner, M., Kellmann, S., Linden, A., Chauhan,S., Funke, B., Grabowski, U., Glatthor, N., Kiefer, M., Schiefer-decker, T., Stiller, G. P., and Versick, S.: Retrieval of temperature, H2O, O3, HNO3, CH4, N2O, ClONO2 and ClO from MIPAS reduced resolution nominal mode limb emission measurements,

Atmos. Meas. Tech., 2, 159175, doi:http://dx.doi.org/10.5194/amt-2-159-2009

Web End =10.5194/amt-2-159-2009 http://dx.doi.org/10.5194/amt-2-159-2009

Web End = , 2009.

Waugh, D. W. and Hall, T. M.: Age of stratospheric air: theory, observations, and models, Rev. Geophys., 40, 1010, doi:http://dx.doi.org/10.1029/2000RG000101

Web End =10.1029/2000RG000101 http://dx.doi.org/10.1029/2000RG000101

Web End = , 2002.

Wofsy, S. C., Boering, K. A., Daube, Jr., B. C., McElroy, M. B., Loewenstein, M., Podolske, J. R., Elkins, J. W., Dutton, G. S., and Fahey, D. W.: Vertical transport rates in the stratosphere in 1993 from observations of CO2, N2O and CH4., Geophys. Res.

Lett., 21, 25712574, 1994.

Wunsch, C.: The ocean circulation inverse problem, CambridgeUniversity Press, 1996.

Atmos. Chem. Phys., 16, 1456314584, 2016 www.atmos-chem-phys.net/16/14563/2016/