Determination of the source regions for surface to stratosphere transport: An Eulerian backtracking approach

A wide range of inverse problems in atmospheric transport and chemistry can be solved within the Eulerian backtracking framework. Here it is shown how a new and accurate numerical implementation can be used as an alternative to Lagrangian back trajectory methods in a wide class of process studies. As a key example, the question of how the (time‐averaged) stratospheric flux of a finite lifetime chemical species depends upon the location(s) of its surface source(s) is addressed. The resulting sensitivity maps are demonstrated to be robust features of the global atmospheric circulation, with relatively low interannual variability. The maps serve as an at‐a‐glance resource for policymakers wishing to compare the likely impact of proposed emission locations for very short lived halogenated species on the total loading of stratospheric chlorine and bromine.


Introduction
A wide range of problems in atmospheric transport and chemistry can be modeled, subject to suitable boundary conditions, by the linear equation Here c ( , t) is the mass mixing ratio of the relevant trace gas, s( , t) is its mass source, ( , t) the density of air, and  is the linear advection-diffusion-reaction-convection operator defined by Here ( , t) is the local mean wind speed, ( , t) a symmetric eddy diffusivity tensor, l( , t) is the local loss rate, e.g., due to photolysis or reaction with a reservoir species, and  is a linear operator modeling the nonlocal transport associated with unresolved convection. Note that (1) can easily be extended to multiple species by replacing c and s by vectors and l by a matrix, which could be, for example, the tangent linear model to the chemistry scheme used by a chemistry transport model (CTM).
A typical objective in solving (1) is to evaluate an integral quantity , e.g., the (weighted) average of c over a given region and time period. The question of interest may then involve determining the sensitivity of  to different configurations of the source distribution s. It is then well known [e.g., Enting, 2002] that rather than solve a large number of forward problems each with different s, it is more efficient to solve the adjoint or inverse equation to (1). Consequently, the theory, development, and numerical implementation of adjoints models for CTMs has been the subject of much research [e.g., Vukicevic and Hess, 2000;Henze et al., 2007].
Many operational adjoint models are derived by (automated) differentiation of the discretized version of the forward model. Eulerian backtracking  provides an alternative based on deriving and then solving the retro transport equation, described in section 2 below, corresponding to (1). The advantages of Eulerian backtracking include 1. The conceptual framework for inverse problems corresponds closely to the "back trajectory" framework exploited in Lagrangian inverse problems [e.g., Seibert and Frank, 2004]. In contrast to other adjoint formulations, Eulerian backtracking inverse problems are therefore simpler to define, understand, and compare with Lagrangian results. 2. The numerical transport scheme used to solve the retro transport equation is essentially the same as that utilized by the forward model. The qualitative behavior of numerical solutions is therefore well understood, and possible numerical stability problems [e.g., Sirkes and Tziperman, 1997] associated with adjoints are avoided.
Hitherto, the principal disadvantage has been that it is far from straightforward to obtain a high level of numerical accuracy in Eulerian backtracking calculations, e.g., , comparing direct and adjoint sensitivities for a short-time test problem, report relative errors of order 10 −2 . To address this, a new and accurate Eulerian-backtracking model (RETRO-TOM) [Haines et al., 2014], based on the transport component of the CTM TOMCAT [Chipperfield, 2006] has been developed and tested by the authors. Relative accuracies of order 10 −8 are reported for RETRO-TOM in comparable test problems.
In section 2, RETRO-TOM is shown to be an efficient alternative to, and a valuable numerical benchmark for, Lagrangian trajectory methods in a wide class of process studies. Problems in which the quantity of interest is either an integral  (as above), or the flux  of a species into a region or air mass, can each be addressed. In section 3, a key example problem addressing the question of how the stratospheric flux  of a finite lifetime chemical species depends on the location(s) of its surface source(s), is revisited using RETRO-TOM. The question is of interest both in quantifying some fundamental climatological transport properties of the general circulation [c.f. Berthet et al., 2007] and of practical value in evaluating the likely contributions of localized sources of very short lived halogenated species (VSLS hereafter) to the stratospheric loading of chlorine and bromine. Source location sensitivities for VSLS have to an extent been previously addressed in forward Eulerian [Wuebbles et al., 2001;Warwick et al., 2006a;Levine et al., 2007;Aschmann et al., 2009;Holzer and Polvani, 2013] and forward Lagrangian [Levine et al., 2007;Brioude et al., 2010;Pisso et al., 2010] studies but not with the flexibility and resolution afforded by Eulerian backtracking. In section 4, conclusions are drawn.

The Eulerian Backtracking Framework
For definiteness, the following quite general form for the convective parameterization  is adopted, following the transilient matrix formalism [Stull, 1984;Romps and Kuang, 2011] Here b(z, z ′ ) represents the nonlocal mass flux density from vertical level z ′ to level z, and is diagnosed from the details of the local meteorological environment at each horizontal position according to a suitable parameterization [e.g., Tiedtke, 1989], and w e is the environmental vertical velocity correction given by A key insight of  is that the Eulerian backtracking framework follows from using the density-weighted inner product to define the adjoint operator  † of . Specifically, if for real-valued functions f and g, and where the spatial integral is over the entire domain Ω, then  † is defined by A straightforward exercise in integration by parts , assuming no-flux conditions at the Earth's surface, reveals that where is the "transpose" of the convection operator .
How can the Eulerian backtracking framework now be used to evaluate the sensitivity of an integral quantity  to changes in the source distribution s? Suppose first that  is the time-integrated total mass of a species HAINES AND ESLER ©2014. The Authors.
in a region  ⊂ Ω, during a time period t 1 < t < t 2 . Then  = ⟨r, c⟩, where the receptor function in this example is If the retro transport equation is then defined to be the definition of the adjoint operator can be used to write The form (11) allows the sensitivity function c * ( , t), which must be found in practice by integrating (10) backward in time, to be interpreted as c * = ∕ s, i.e., proportional to the rate of change of  with respect to a change in the total source s at position and time ( , t). Since the problem (1) is linear, knowledge of c * in the source region is sufficient to obtain  for any source distribution s, simply by evaluating the integral defined by (11).
In a wide range of problems it is the mass flux  of the trace gas into a region  (rather than ) that is the quantity of interest. This case can be addressed by first writing where  is the (time-dependent) bounding surface of , its unit normal, and is the velocity of the moving surface itself (parallel to ). Focusing here only on the first (advective) term in (12), careful application of the divergence theorem allows  to be expressed as an inner product  = ⟨r, c⟩, where in this case In fact, it is established in Text S1 in the supporting information that (13) holds exactly even when the convective and diffusive components of  are treated explicitly. Evidently c * =  ∕ s can now be found by integrating the retro-transport equation (10) forced by the receptor field (13).

Sensitivity Maps of Stratospheric Fluxes of VSLS-Like Tracers
Next, RETRO-TOM is used to calculate the sensitivity of the mass flux  into the stratospheric overworld of a finite lifetime chemical species, to the location of its surface source, which for simplicity is taken to be uniform in time during a given period (t a ≤ t ≤ t b ). In each of the forward problems corresponding to the Eulerian backtracking calculation, the source is located at an isolated point on the Earth's surface and emits a total mass  in during its active period. While a forward integration of (1) could be used to calculate the ratio  ∕ in at a single point, a single integration of (10) produces a "sensitivity map" of  ∕ in at all longitudes and latitudes ( , ) on the Earth's surface, by evaluating the integral where z s ( , ) is the surface altitude.
The receptor region  in (13) is taken to be the stratospheric overworld (potential temperature > 380 K), following Pisso et al. [2010] who have demonstrated its relevance for calculating ozone depletion potentials (ODPs) of VSLS. If desired, a latitudinal weighting on the flux  , to account for the variations in stratospheric residence times of the active halogen species noted by Pisso et al. [2010], could easily be introduced. However, Pisso et al. [2010] observed that it is primarily the stratospheric residence time associated with tropical injection that determines ODP and so, in the interests of simplicity, we keep a uniform weighting here. The receptor is active once the emissions commence (t 1 = t a , t 2 → ∞), and the chemical tracer itself is chosen to have a uniform lifetime l( , t) = 1∕ . Again, this is an idealization; however, Brioude et al. reported for several VSLS that  is considerably more sensitive to emission location than to spatial variations in chemical lifetime.
Equation (10) where T ≥ 1 year, is sufficient to get good convergence of results). The convective transport operator b(z, z ′ ) is obtained using the convective parameterization of Tiedtke [1989] and explicit diffusion is nonzero only in the boundary layer where the parameterization scheme of Louis [1979] is used. A discussion of the biases introduced by these schemes is given by Hoyle et al. [2011]. Full details of the RETRO-TOM implementation are provided by Haines et al. [2014]. Figure 1 shows sensitivity maps calculated for a species with = 20 days emitted during a 3 month season. Results are averaged over 4 years (December 2007to November 2011. In each case surface to stratosphere transport of the species is clearly dominated by emissions from the tropics (average value of  ∕ in = 1%) as opposed to the extratropics ( ∕ in = 0.2%). These averaged values are about ten times lower than those reported in forward calculations by Holzer and Polvani [2013], but the discrepancy can be explained by our use of the stratospheric overworld rather than the entire stratosphere as the receptor. Consistent with this finding Levine et al. [2007] report that, of trajectories released from the surface subsequently entering the stratosphere, approximately 95% enter the lowermost stratosphere (bounded below by the 2 potential vorticity unit surface and above by = 380 K) with only 5% entering the overworld ( > 380 K).
Our results have been validated by comparison with a set of 16 forward integrations of (1) covering the period June 2011 to June 2012. In each forward run, the source is confined to a single grid-box on the 138 • longitude circle, with latitudes 7 • S-35 • N (see Figure 2c, red line). The quantity  ∕ in is calculated explicitly for each forward run using TOMCAT with the results shown in Figure 2a. Excellent agreement is obtained, with relative errors of order 10 −8 . Note that in the forward runs flux limiters are turned off in the advection scheme, which both leads to more accurate forward advection (as measured by a global norm) and is required for close agreement with RETRO-TOM. However, such an approach is not suitable for all chemistry applications as small regions of negative mixing ratio can result (see discussion in Haines et al. [2014]). The strong seasonal variation in  ∕ in shown in Figure 1, notably its maxima over the Himalayan Plateau and the Indian subcontinent in JJA and over the maritime continent in DJF, has been noted previously [e.g., Berthet et al., 2007;Pisso et al., 2010], although not with the detail shown here. A key question, easily answered using RETRO-TOM, concerns the extent to which such detailed pictures are meaningful, i.e., how much do the sensitivity fields in Figure 1 vary from year to year? Figures 2b and  2c show the sensitivity maps for constant emission of a 20 day lifetime tracer in JJA 2010 and 2011, respectively. The remarkably low interannual variability seen in Figures 2b and 2c is repeated across all seasons and years. The sensitivity maps in Figure 1 are therefore robust features of the tropical climatology, insensitive to weather fluctuations, and with at most weak dependence upon global climate indices such as El Niño-Southern Oscillation (the period of the study covers the 2009-2010 El Niño and several distinct La Niña events), with the caveat that other effects (e.g., volcanic eruptions) might be revealed by a longer study.
We have performed several calculations to confirm that the above findings are not sensitive to details of the model formulation. First, the JJA 2011 sensitivity map calculation (shown in Figure 2c) was repeated at the increased horizontal resolution of 1.125 • × 1.125 • , with the results shown in Figure 2d. Other than a slight increase in the average flux (22%) the sensitivity pattern is almost unchanged. Second, we have repeated our calculations with the convection parameterization switched off (not shown). For species with ≳ 20 days, differences in  ∕ in are surprisingly small, with the largest changes being increased transport from the western Pacific region in JJA (approximately 40% greater) and decreased transport over the Himalaya but with the overall pattern relatively unchanged. Recent studies [Russo et al., 2011;Hoyle et al., 2011], however, indicate that the Tiedtke scheme underestimates cloud top heights in convectively active regions such as the western Pacific, and consequently, our "convection on" experiments may be underestimating transport there.
To illustrate how sensitivity maps vary with tracer lifetime, Figure 3 shows results for tracers with = 10, 20, 50, and 100 days, corresponding to forward runs for which the source is active for all 4 years of our study. The global average of  ∕ in increases from 0.16% for the =10 days tracer to 2.7% for =100 days. The corresponding figures for the tropical average of  ∕ in are 0.4% and 3.8%, and for the extratropics, 0.05% and 2%. Clearly, as tracer lifetime increases, transport to the stratosphere from extratropical sources HAINES AND ESLER ©2014. The Authors. becomes proportionately much more significant. Zonal variations in  ∕ in are also seen to become much weaker as increases, and for the = 100 day tracer  ∕ in is well approximated by a function of latitude only. By contrast,  ∕ in for the = 10 day tracer is large only in the western Pacific, Indian subcontinent, and to a lesser extent in equatorial South America.

Conclusions
It has been demonstrated above that the Eulerian backtracking method of  can be used as a flexible and efficient alternative to Lagrangian back trajectories in numerous problems of interest to atmospheric scientists, in particular those where the quantity of interest can be expressed as an integral  or a flux  . For example, further uses of RETRO-TOM could include evaluating source regions for pollutant transport into the Arctic boundary layer in winter [Stohl, 2006], or establishing the extent to which upstream emissions at different possible sites influence localized aircraft measurements far downstream [e.g., Stohl et al., 2003]. The results obtained would serve to corroborate and extend the Lagrangian studies cited, within a deterministic framework that is free from statistical error associated with stochastic trajectories. The potential for coupling RETRO-TOM with the adjoint to a tangent linear model of an atmospheric chemistry scheme is also clear.
Here the key example problem of transport from the surface to the stratosphere of VSLS-like tracers has been investigated. The results from RETRO-TOM have been expressed as sensitivity maps of  ∕ in , the ratio of emitted tracer reaching the stratospheric overworld from every possible surface emission site. RETRO-TOM has revealed that these maps, which are subject to a strong seasonal cycle, are subject to remarkably little interannual variability. The routes for surface to stratosphere transport, at least on the seasonal time scale, are found to be robust features of the (tropical) circulation. Figure 3 can therefore serve as an at-a-glance guide to the proportion of emitted chlorine reaching the stratosphere for different potential new anthropogenic emission sites of (say) a chlorinated VSLS, which may be of interest to policymakers. Species-specific sensitivity maps, with a more sophisticated treatment of VSLS chemistry, can be generated with just a few modifications to RETRO-TOM.