Interactive comment on “ Technical Note : Adjoint formulation of the TOMCAT atmospheric transport scheme in the Eulerian backtracking framework ( RETRO-TOM ) ” by P . E

2. We have not attempted a test case with increased vertical resolution, as the vertical resolution is significantly more difficult to change in TOMCAT. However, we do not see any reason why increasing the vertical resolution would cause any loss of accuracy of our algorithm. The column-matrix formulation of the parametrisations means that no essential change to the adjoint formulation is necessary when vertical resolution is increased.


Introduction
The past 20 years or so have seen an explosion in the development of adjoint models for chemistry transport models (CTMs).Adjoint models have numerous applications (e.g.Enting, 2002), including variational data assimilation of constituent concentrations (Elbern and Schmidt, 1999), inverse modelling of chemical source strengths (Müller and Stavrakou, 2005;Meirink et al., 2006), sensitivity analysis (Vukicevic and Hess, 2000) and parameter sensitivity esti-mation (e.g.Navon, 1998).The numerical accuracy (defined here as the relative difference between sensitivities calculated by a linear perturbation to a forward calculation and those obtained from the adjoint) as well as the reliability of adjoint models are evidently key to the above applications.Data assimilation applications in particular are sensitive to numerical inconsistencies between the formulations of the forward and adjoint models.Careful numerical analysis in the development of adjoint models is therefore crucial.
The focus of this work will be the linear advectiondiffusion-convection problem that forms the "dynamical core" of a CTM, although we believe that all of our results can be straightforwardly extended to the case of the tangent linear model of the full CTM.For the related problem of nonlinear general circulation models, two approaches to formulating adjoints are summarized by Sirkes and Tziperman (1997).The methods are namely "finite difference of adjoint" (FDA) and "adjoint of finite difference" (AFD).In the context of the linear problem considered here, the FDA approach involves finding the (continuous) adjoint equation for the underlying (continuous) advection-diffusion equation, followed by discretizing the continuous equation.In the AFD approach, the forward equation is first discretized, and then the adjoint of the resulting discrete system of equations is taken.
The advantage of the FDA approach is that a partial differential equation is obtained, which can then be solved by reliable and well-understood numerical methods.A disadvantage is possible numerical inaccuracy, of the order of the discretization error, compared with sensitivities calculated using the forward model.The advantage of the AFD approach Published by Copernicus Publications on behalf of the European Geosciences Union.
is that the instantaneous sensitivities calculated by the adjoint model typically match those of the forward model to within machine precision.The disadvantage is that, as emphasized by Sirkes and Tziperman (1997), the time-stepping behaviour of the adjoint equations is poorly understood, and it is possible that spurious computational modes overwhelm the calculations.In particular, several authors (e.g.Hourdin et al., 2006;Henze et al., 2007;Gou and Sandu, 2011) have noted that AFD adjoints to nonlinear advection schemes can lead to undesirable results and poor performance.
The aim of the present work is to describe the development of an adjoint RETRO-TOM to the dynamical core of the TOMCAT CTM (Chipperfield, 2006) that combines the desirable numerical and conceptual properties of the FDA approach with the accuracy of an AFD model.Our results are achieved by utilizing the Eulerian backtracking framework of Hourdin and Talagrand (2006), which is an elegant variation on the FDA approach which maximizes the symmetry between forward and adjoint models.Specifically Eulerian backtracking confers the following advantages: -There is a close correspondence between calculations made using the retro-transport equation of the Eulerian backtracking method and Lagrangian back trajectory calculations (e.g.Seibert and Frank, 2004).Compared with alternative formulations, Eulerian backtracking results are simpler to define, understand and compare with Lagrangian results, allowing for their ease of use in process studies.
-The numerical transport scheme used by the Eulerian backtracking model is essentially that utilized by the forward model.The qualitative behaviour of numerical solutions is therefore well understood, and numerical stability problems such as those discussed above can be avoided.
One reason to suspect a priori that an accurate adjoint of TOMCAT could be formulated using Eulerian backtracking is that the Prather (1986) advection scheme used by TOM-CAT can be shown (under certain circumstances, see below) to be time symmetric in the sense of Hourdin et al. (2006).In other words, the Prather scheme applied to the retro-transport equation turns out to be the exact numerical adjoint of the Prather scheme applied to the forward problem.Alternatively, in the language of Sirkes and Tziperman (1997), the finite difference of the adjoint coincides with the adjoint of the finite difference so that the desirable numerical properties of the FDA are combined with the accuracy of the AFD.The column-matrix formulation of both the convective and boundary-layer turbulence parameterizations in TOMCAT also lend themselves to accurate adjoint formulation.
Notwithstanding the above, the key development that enables the high accuracy obtained by RETRO-TOM is a careful treatment of what might be termed the "density incon-sistency problem" of offline forward CTMs.First discussed in detail by Jöckel et al. (2001), density inconsistency arises when forcing wind and density fields (obtained in TOMCAT from surface pressure) are provided (e.g. from re-analysis products) at finite time intervals.There is then an inconsistency between the density field computed by forward advection from the previous forcing field and that obtained from the new forcing field.TOMCAT addresses density inconsistency by making a discontinuous update to the density field as each new forcing file is read.In order that chemical species in TOMCAT remain "well-mixed" (in the sense that spatially uniform mixing ratios remain so), the density update is also applied to tracer mass fields meaning that global mass conservation of each tracer is violated.Because the issue is fundamental to the intermittent nature of the forcing files, other CTMs must unavoidably address density inconsistency in a similar way.
Our position is that a CTM adjoint should be built upon a numerical scheme for the "dynamical core" that is both highly accurate and numerically well formulated, at least in the absence of advection-scheme-related nonlinearities due to e.g.flux limiters (Thuburn and Haine, 2001;Vukicevic et al., 2001;Hourdin et al., 2006) which raise various separate issues.Such a model provides as solid as possible a foundation for the applications listed above, and in particular allows numerical errors to be excluded, as far as possible, when assessing results.Previous efforts to assess the accuracy of the transport component of CTM adjoint calculations have been made by Hourdin et al. (2006), Henze et al. (2007) and Wilson et al. (2013).Hourdin et al. (2006, see their Fig. 2b) record relative errors, comparing adjoint sensitivities and sensitivities calculated directly using the forward model in a 3-day integration, of the order of 0.5 % when the linear Godunov (upwind) scheme is used.Henze et al. (2007, see their Fig. 7) show significantly larger errors for a 2-day integration, although it is likely (see referee comment by D. Henze) that their poorer results are predominantly due to the difficulties of generating an adjoint to the nonlinear piecewise parabolic scheme (Lin and Rood, 1996) used for their test case.Both studies suggest significant room for improvement that has motivated the present work, which develops an adjoint model using Prather's advection scheme.The recent AFD scheme of Wilson et al. (2013) also uses Prather's scheme, and has demonstrated accuracies close to machine precision, but could be subject to the problems of AFD adjoints to nonlinear advection schemes highlighted by e.g.Gou and Sandu (2011).The present study is the first to exploit the time-symmetric properties of the Prather scheme.
Eulerian backtracking, i.e. the framework of Hourdin and Talagrand (2006) underpinning RETRO-TOM, is summarized in Sect. 2. In Sect.3, the formulation of RETRO-TOM is described in detail.Section 3.1 details the key aspects of the forward model (TOMCAT) including the treatment of density inconsistency.Section 3.2 continues by demonstrating the time symmetry of the Prather (1986) scheme and Sect.3.3 gives the explicit numerical details of RETRO-TOM.In Sect. 4 the accuracy of RETRO-TOM is assessed in three different test problems, in which the transport characteristics are dominated by short-term advection (Sect.4.1), short-term convection (Sect.4.2) and long-term inter-hemispheric transport (Sect.4.3), respectively.The extent to which the accuracy of RETRO-TOM depends upon a correct treatment of the "density inconsistency" problem is discussed in Sect.4.4.The difficulties associated with flux limiters are discussed in detail in Sect.4.5, and the effect of numerical resolution of the forward model is considered in Sect.4.6.Conclusions are drawn in Sect. 5.

Mathematical formulation of Eulerian backtracking
We consider a transport problem described by a linear equation of the form (1) Here c(x, t) is the mixing ratio of the relevant trace gas, s(x, t) is its source, ρ(x, t) the density of air and L is the linear advection-diffusion-reaction-convection operator defined by Here u(x, t) is the local mean wind speed, κ(x, t) is a symmetric eddy diffusivity tensor, l(x, t) is the local loss rate e.g.due to photolysis or reaction with a reservoir species, and C is a linear operator modelling the non-local transport associated with unresolved convection (see Haines and Esler, 2014, for an exact definition of C).
The Eulerian backtracking formulation (Hourdin and Talagrand, 2006) follows from using the density-weighted inner product to define the adjoint operator for real-valued functions f and g, and where the spatial integral is over the entire domain , then L † is defined by A straightforward exercise in integration by parts (Hourdin and Talagrand, 2006), assuming no-flux conditions at the Earth's surface, reveals that where C † is the "transpose" of the convection operator C.
A key insight of Hourdin and Talagrand (2006) is that using the density-weighted inner product (3) to define the adjoint operator leads to a symmetry between the direct and adjoint equations.As a result of this symmetry the adjoint, or retro-transport, equation corresponds to solving the forward problem backwards in time and with the advective mass fluxes reversed.Other researchers (e.g.Sandu et al., 2005;Hakami et al., 2007;Henze et al., 2007;Gou and Sandu, 2011) have used non-density-weighted inner products to construct a continuous adjoint model which can then be discretized and solved.The result is an asymmetry between the form of the forward transport equation (e.g. an advection equation written in terms of mixing ratio) and its adjoint (e.g. a flux-form conservation law written in terms of tracer mass).Compare, for example, Eqs. (1a) and (5a) of Sandu et al. (2005).The disadvantage, apart from inelegance, of this approach, is that if the same numerical scheme is used to solve both equations (Henze et al., 2007), then numerical inaccuracies can be introduced in moving between one form and the other.For example, Henze et al. (2007) rescale their adjoint variable before and after each advective time step, by dividing and multiplying by the air density respectively, under the approximation that the density remains constant across each time step.See also Table 2 of Hakami et al. (2007) which details the steps required to convert their adjoint variable to and from a mixing ratio during each model step.
In practice, the ultimate object of solving Eq. ( 1) is often to evaluate an integral quantity I = r, c with r(x, t) a "receptor" function defining the location and time of the measurement, and allowing for the possibility of non-uniform spatial and temporal weighting.In many process studies, the question of interest involves determining the sensitivity of I 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 one adjoint or inverse equation to (1) from which the same sensitivity of I to s can be evaluated.If the retro-transport equation is defined to be where c * is the mixing ratio of a "retro-tracer", then the definition of the Eulerian backtracking operator L † can be used to write The form (7) allows I to be calculated from the retro-tracer c * (x, t) obtained by solving Eq. ( 6), which must be found in practice by integrating backwards in time.The retrotracer c * (x, t) is, equivalently, the sensitivity of I with respect to a change in s at the location and time (x, t), as can be expressed mathematically by the functional derivative c * = δI/δs.Since the problem (1) is linear, knowledge of c * throughout the source region is sufficient to obtain I for any given source distribution s, simply by evaluating the integral defined by Eq. ( 7).

Numerical formulation of TOMCAT's forward scheme
TOMCAT is a global three-dimensional off-line chemistry transport model, which is run here with Gaussian horizontal grids of size 128 × 64 (approx.2.8 • × 2.8 • ) and 320 × 160 (approx.1.1 • × 1.1 • ).A hybrid sigma-pressure coordinate is used in the vertical direction with 31 model levels extending from the surface up to approximately 30 km.Advection is performed using the second-order moments scheme of Prather (1986) with the forcing wind fields obtained from the European Centre for Medium-Range Weather Forecasts (ECMWF) operational analyses.TOMCAT parameterizes subgrid-scale deep convection according to the scheme of Tiedtke (1989) and boundary layer vertical diffusion according to Louis (1979) (see Stockwell and Chipperfield, 1999, for further details).An assessment of the performance of the convective parameterization in TOMCAT is given by Hoyle et al. (2011).Typically a new forcing file is supplied once every 6 h (the "forcing period" hereafter) and the required mass fluxes, temperature field etc. are linearly interpolated across the forcing period.
A brief outline of TOMCAT's dynamical core, as required to understand the formulation of RETRO-TOM, is described next for a single species subject to chemical loss.The rectangular grid in longitude and latitude and the 31 vertical levels together divide the atmosphere into a total of N grid cells.It is helpful to express quantities such as the tracer mixing ratio and the air mass in each box as vectors of length N, i.e. c = (c 1 , . . ., c N ) T and m = (m 1 , . . ., m N ) T .Further, the total mass of tracer in each grid box can be written S = (S 1 , . . ., S N ) T  = Mc, where M = diag(m) is an N × N diagonal matrix containing the box masses.
The discrete form of the transport operator Eq. (1) in TOMCAT, at the nth time level, can be written as an operator L n satisfying c n+1 = L n c n .The discretized operator L n can be further decomposed into a successive application of a number of sub-operators representing advection, chemical loss and non-local vertical transport.The advection sub-operator can be further split into individual components that perform advection in each coordinate direction.Finally, TOMCAT's treatment of the density inconsistency problem requires a density correction operator to be applied at the start of each forcing period, In the Prather scheme, the advection operator in fact acts on a longer state vector of length 10N, consisting of the total tracer mass and first and second moments of the tracer in each grid cell.For simplicity we describe the operators here in terms of their action upon the first of these 10 components, the zeroth-order moment giving the total mass of tracer S in each grid cell.With the exception of advection, which will be considered in detail in Sect.3.2 and in Appendix A, higherorder moments are treated in the same way as the zeroth-order moment.We now summarize these operators, adopting the convention that − and + subscripts indicate the value of a variable respectively before and after the application of a particular operator.
Subgrid-scale convection and boundary layer vertical diffusion are implemented by multiplying each vertical column of the grid by an individual matrix operator.For mathematical convenience, we represent this process in terms of a single operator V n acting on all the tracer masses S simultaneously, For computational efficiency TOMCAT calculates the operator V n only once per forcing period.In order to ensure that a tracer distribution corresponding to a uniform mixing ratio remains so under the action of V n , it is necessary to also apply the operator V n to the box masses m as part of each convective step: An operator D n accounts for chemical loss over each time step ( t), where l i,n = l(x i , t n ) gives the local loss rate in the grid cell centred on x i at time t n .
As described in the introduction there is an inconsistency between the density field obtained by advection over a forcing period and the density field implied by the surface pressure supplied by the next forcing file.TOMCAT updates the density field at the start of each forcing period so that it is in agreement with the surface pressure supplied by the new forcing file.In order to preserve the mixing ratio c + = c − across this update the tracer mass in each grid box is adjusted simultaneously.An instantaneous update to box masses from m − to m + requires that we update the tracer masses according to at the beginning of each forcing period.

Time symmetry of the Prather scheme
In this section we discuss the time symmetry of the Prather (1986) advection scheme used by TOMCAT.A numerical scheme is described as time symmetric if the matrix operator generated by its application to the retro-transport equation is an exact adjoint of that generated by its application to the forward problem (Hourdin et al., 2006).Time symmetry has been previously demonstrated in general for the first-order upwind Godunov scheme (Hourdin et al., 2006), and for the quadratic upstream interpolation algorithm (QUICK) in the special case of one-dimensional advection by a non-divergent (i.e.spatially uniform) wind field (Vukicevic et al., 2001).The Godunov scheme is equivalent to advection by zeroth-order moments.
Our discovery that the Prather scheme is time symmetric in general motivated the development of RETRO-TOM because of the promise of high accuracy.Here and in Appendix A an explicit proof is presented for the time symmetry of the first-order moment scheme (Russell and Lerner, 1981), which is equivalent to the Prather scheme with second-order moments neglected.The full exposition of the proof for the Prather scheme itself is too lengthy to be given in full here, but proceeds by exact analogy with the proof for the firstorder scheme (however, one key result for the Prather scheme is given in Appendix A5 by means of illustration).An important detail is that it is the mass flux (rather than the wind speed) that is reversed in the adjoint calculation.
Advection in higher dimensions is implemented in the Prather scheme by successively applying the operators associated with one-dimensional advection in each coordinate direction (time splitting, Strang, 1968), hence time symmetry need be established for one-dimensional advection only.The time symmetry of higher-dimensional advection follows by treating the coefficients associated with variation in the other dimensions in terms of their variation in the single dimension in which advection is taking place.
As discussed in Sect.3.1 the Prather scheme acts upon a total of 10 components per grid cell.Introducing an extended state vector s containing these moments, defined precisely in Appendix A, the advection operator at time t n can be expressed as s + = A n s − .Here the − and + subscripts indicate variables evaluated before and after application of A n respectively.Following Hourdin et al. (2006) we define a discretized inner product, evaluated at a single time, by where f and g are column vectors and the "weighting matrix" W is a diagonal matrix constructed from the box air masses and described in Appendix A. The adjoint of an operator A with respect to Eq. ( 12) must then satisfy A † f , g − = f , Ag + for all possible f and g.That is Here we drop the subscript n to indicate that the result (13) applies equally well to both the umbrella advection operator A n and to the individual advection operators in each dimension, provided of course that the weighting matrices W − and W + are constructed at the appropriate point.
The RETRO-TOM advection operator B n for step n is obtained by discretizing the advective part of Eq. ( 6) using the same Prather scheme as in the forward model.The operator B n is applied to the extended state vector s * associated with Eulerian backtracking according to s * − = B n s * + .A numerical scheme is time symmetric if the forward operator A and Eulerian backtracking operator B satisfy In Appendix A we construct the matrix operator A for one-dimensional advection by first-order moments and then demonstrate explicitly that A and B satisfy time symmetry (i.e.Eq. 14).A brief outline of the second-order moment result is then provided.
A simple summary of the method follows.The major part of the work is in combining the "splitting" and "recombining" stages of Prather's algorithm into a single matrix operator, from which the components of A and B can be obtained.We start by describing the basis functions and coefficients associated with approximating a one-dimensional function by first-order moments, and then demonstrate how they are manipulated.These results are then applied to advective transport by firstly splitting each grid box into three parts.The first part contains the tracer to be advected to the left, the second part remains in the original grid box, and the third part is advected to the right.These sub-boxes are then combined into new boxes, each with a single set of moment coefficients.The equations for splitting and recombining moments are then combined into a single set of expressions giving new moments for grid box i in terms of the old moments in grid boxes i − 1, i and i + 1 in a single step.Finally, the equations obtained are converted into a matrix stencil for which Eq. ( 14) can be verified.

Numerical implementation of RETRO-TOM
The time symmetry of the Prather scheme is conditional on the density field for Eulerian backtracking being identical to that used in the forward run.Due to the "density inconsistency" between surface pressure and forcing wind fields discussed in Sect. 1 the density field obtained by TOMCAT after a forward integration over a forcing period will differ from that supplied by the forcing file associated with the next forcing time.If this second density field is used as the starting point for the (backwards) RETRO-TOM calculation across the forcing period there will be an associated discrepancy between the forward and Eulerian backtracking density fields.To avoid this discrepancy the density field in RETRO-TOM is instead initialized at the start of each forcing period to agree with that obtained in the forward problem at the same time.
To initialize the density field at the start of each forcing period in RETRO-TOM, the box masses are first advected forwards across the forcing period, exactly as in a TOM-CAT calculation, in order to obtain the appropriate density field.This additional forward calculation does not lead to a significant increase in the length of the run as the forward transport of the density field is an order of magnitude lower in computational cost compared to a single tracer species (since there are 10-second-order moments per species).The resulting density field is stored at every model time step between the two forcing times, for use with the RETRO-TOM backwards integration.The above procedure ensures that RETRO-TOM uses an identical density field to TOM-CAT at every model time step.
While the advective adjoint operator is obtained from a discretization of the retro-transport Eq. (A23), exact adjoints to the other operators that comprise TOMCAT are more readily obtained directly from their forward counterparts.As seen in Sect.3.1, individual operators in TOMCAT and RETRO-TOM are best described in terms of their effect upon the total mass of tracer in each grid cell, S. The discrete inner product associated with Eq. ( 3) is then where F and G are column vectors of length N. The inner product Eq.( 15) is used to define adjoints to the forward transport operators.For example, the adjoint of the discretized convective operator V n (satisfying Eq. 8), is given by In RETRO-TOM, V † n acts upon the Eulerian backtracking equivalent of S, the total mass of retro-tracer in each grid box S * , according to S * − = V † n S * + .The forward convective scheme also changes the density field by a small amount (see Eq. 9) to ensure that the convective scheme preserves a uniform mixing ratio.Although this density update is small in comparison with that associated with reading a new forcing file, in order to enable RETRO-TOM to proceed with an identical density field to the forward problem it is necessary to reverse this update.The relevant box masses are already calculated when advecting the box masses forward across each forcing period.Thus, although proscribing the density at each time step introduces a slight storage overhead, it does not represent a significant computational one.
The adjoint to the operator Eq. ( 11), which preserves a uniform mixing ratio across an instantaneous change in air density, is the identity matrix That is, no adjustment to S * is required in RETRO-TOM in response to an instantaneous change in air density.
Since the tracer decay operator D n is diagonal and does not alter the mass matrix M it remains as for the forward problem.To further increase the accuracy of the adjoint model the order in which the individual operators are applied is reversed in RETRO-TOM.At each time step the order advection-decay-vertical mixing in TOMCAT becomes vertical mixing-decay-advection in RETRO-TOM.It will be demonstrated below that, with the correct density field proscribed and the adjoint operators implemented as described here, RETRO-TOM yields an accurate adjoint to TOMCAT.

Validation of RETRO-TOM
In order to validate RETRO-TOM, three separate model problems are considered next.These problems are intended to examine different aspects of the model functionality, namely (I) short-term primarily advective transport in the extratropics, (II) short-term primarily convective transport in the tropics, and (III) inter-seasonal inter-hemispheric transport.RETRO-TOM will be validated by a comparison between sensitivities calculated in a single Eulerian backtracking calculation, and sensitivities obtained by a number of forward calculations with isolated sources.
In each case we are interested in the sensitivity of a timeintegrated quantity of tracer (I) over a region D of the atmosphere, with respect to surface emissions at some earlier time.In the inner product notation introduced in Sect.2, I = r, c is the integral, with respect to space and time, of tracer mass ρc over the region D and the detection period (t 1 , t 2 ).This corresponds to the receptor function r(x, t) being For definiteness, we aim to calculate the sensitivity of I to the surface emissions during an earlier period (t a < t < t b ).
The result can be expressed as a map for longitude λ, latitude φ, surface altitude z s (λ, φ) and with c * = δI/δs obtained by using RETRO-TOM to solve Eq. (6) with r(x, t) as in Eq. ( 18).
Alternatively, the map C * (λ, φ) can be constructed "grid cell by grid cell", by solving Eq. ( 1), in a (possibly large) number of forward TOMCAT calculations.In each of these calculations the tracer source s(x, t) is confined to a single grid cell and I is calculated by integrating c(x, t) over the receptor region (i.e.calculating the inner product r, c directly for each grid-cell-sized source.The forward sensitivity map is built up systematically from these integrations, a method that is obviously highly inefficient compared with the single RETRO-TOM calculation described above.The results presented below are for TOMCAT and RETRO-TOM run at a horizontal resolution of 128 × 64; a comparison with results obtained at higher resolution (320×160) is given in Sect.4.6.

Test case I
Test case I is designed to test the accuracy of RETRO-TOM's advective transport operator, by solving a relatively short-term primarily advective transport problem in the extratropics.The time period of interest is 1-11 January 2012 (00:00-00:00 UT).The receptor region D is centred over central Europe (13 • < λ < 24 • and 42 • < φ < 53 • ) and situated in the upper troposphere at a height of approximately 6.5-9 km.The receptor region D corresponds to 4 × 4 × 4 TOMCAT grid cells.Both the emission period and the detection period correspond to the full 10-day run.
Figure 1a shows the sensitivity map C * obtained from a single run of RETRO-TOM.The map shows that the receptor region D (white rectangle in Fig. 1) is most heavily influenced by surface sources in the western Atlantic compared with other locations.To validate the results, Fig. 1a is compared with results obtained by repeatedly solving the forward problem, focusing on a subset of the western Atlantic region demarcated by the red rectangle (10 × 10 surface grid cells).Evidently this requires 100 forward TOMCAT runs.The results of these are shown in Fig. 1b and are indistinguishable by eye from those within the red rectangle in Fig. 1a.
To quantify the difference between the two calculations of C * , we define the local relative error where the subscripts D and R refer to "direct" and "retro" respectively.Figure 1c shows E(λ, φ) with a logarithmic colour scale.It is evident that the sensitivities obtained by Eulerian backtracking have maximum relative error of order 10 −8 .
Figure 1d shows the relative error for a comparison between the direct sensitivities and those computed by Eulerian backtracking without reversing the sequence of individual operators in RETRO-TOM.The maximum and average relative errors over the 10 × 10 patch of grid cells shown are 0.014 and 0.004 respectively.For their 4-day test case Hourdin et al. (2006) reported a decrease in error of the same order when they reversed the sequence of individual operators.Figure 1e will be discussed below.
To assess the importance of convection in test case I, it was repeated with the convective parameterization (although not the boundary layer diffusion) switched off.The results (not shown) showed a modest reduction in the reported sensitivities of approximately 10 %.

Test case II
Test case II is designed to test the accuracy of RETRO-TOM's treatment of convective transport, by solving a relatively short-term vertical transport problem in the tropics.
Here the time period of interest is the first 20 days of January 2012, with the receptor active throughout, and the receptor region D (10 × 10× 4 model grid cells) located over the maritime continent (97 • < λ < 125 • , −14 • < φ < 14 • ) at approximate height 11-14 km, i.e. in the tropical tropopause layer.Unlike in test case I, the source is active for only the first day of the calculation (1 January).
Figure 2a shows the sensitivity map C * to emissions during the 1 January.The map shows that the receptor region D (again indicated by a white rectangle) is principally influenced by surface sources almost directly below, strongly suggesting the importance of deep convection in driving the transport.The region for the forward runs is chosen accordingly (red rectangle).Again 100 forward calculations are required to obtain C * D as shown in Fig. 2b with the direct sensitivities again indistinguishable from those in Fig. 2a. Figure 2c shows the relative error, E, of the sensitivities obtained by Eulerian backtracking with respect to those obtained by direct transport.The maximum relative error is again found to be of order 10 −8 .Fig. 2d-e will be discussed below.
Test case II was repeated with the convective parameterization switched off, and this time the results were profoundly different, with sensitivities approximately 60 % lower than those reported above.

Test case III
The magnitude of the errors in test cases I and II are sufficiently small to motivate a much longer integration to provide a sterner test of the accuracy of RETRO-TOM.Test case III examines inter-hemispheric surface-to-surface transport over a time period of six months (00:00 UT 1 January 2012-00:00 UT 1 July 2012) with both the sources and the receptor active throughout.grid cells centred over the United States.In order to evaluate the ability of RETRO-TOM to correctly incorporate tracer decay the tracer is assumed to decay exponentially with an e-folding time of 50 days.Figure 3a shows values of C * as obtained from RETRO-TOM.The map shows that the influence of Southern Hemisphere sources upon a surface receptor in the Northern Hemisphere is primarily a function of latitude, decreasing by an order of magnitude from north to south between the equator and the pole.Note that a logarithmic colour scale has been used to visualize this large variation.
To assess the accuracy of RETRO-TOM in determining sensitivities for inter-hemisphere transport Fig. 3b compares the sensitivity given in panel a against direct sensitivities obtained for a surface source region made up of a patch of 11×1×1 model grid cells.The patch lies on a latitude circle (38 • S) in the region of New Zealand and the south Pacific and is indicated by the red line on Fig. 3a.Despite the significantly increased timescale the maximum relative error in the value of C * obtained from RETRO-TOM (black line in Fig. 3b in comparison to the direct sensitivity (blue triangles) is again found to be of order 10 −8 .

Importance of the treatment of the "density inconsistency" problem
Here, the extent to which the accuracy of RETRO-TOM depends upon a careful treatment of the "density inconsistency" problem, discussed in Sect. 1, is demonstrated.As stated in Sect.3.3, backwards calculations in RETRO-TOM use the density field calculated in the forward run at each time step.
Here the results obtained are compared with calculations, referred to as "naive RETRO-TOM", in which the densities in RETRO-TOM are updated between forcing files in a similar fashion to the forward model.To be precise, in "naive RETRO-TOM", the forward integration step described in Sect.3.3 is omitted, and the density field is calculated at time steps in between forcing files by solving the mass continuity equation, using the backwards winds.
To quantify the difference between RETRO-TOM and "naive RETRO-TOM" we re-run test cases I and II with the naive formulation.For test case I naive RETRO-TOM results in the errors shown in Fig. 2e.The maximum and average relative errors are 0.005 and 0.001 respectively for the patch shown in Fig. 1a, compared with 10 −8 and 10 −8 in RETRO-TOM.Similarly, for test case II the relative errors are shown in Fig. 2e, with maximum and average values for the patch shown of 0.003 and 0.001 respectively.
A further test involved correctly initializing the density field in RETRO-TOM, following the procedure of Sect.3.3, but not then updating the density field at each time step using the density field calculated from the forward run.In this case errors are introduced by the failure of RETRO-TOM to reverse the effect of the convective operator Eq. ( 9) on the density field.However, the errors in this case are much smaller: a maximum relative error of 1 × 10 −4 is obtained for both test cases I and II.For test case II this error is shown in Fig. 2d.

Evaluation of the effect of flux limiters
The transport problems considered in this work are linear.However, chemistry schemes require positive concentrations and the means by which this is achieved in the Prather scheme is by introducing flux limiters (specifically that described in Sect. 4 of Prather, 1986), which are necessarily nonlinear in tracer concentration.Other widely used schemes (e.g.Lin and Rood, 1996;Hourdin and Armengaud, 1999) use similar devices.It is highly questionable whether it is even desirable to generate an exact adjoint to a nonlinear advection scheme.Several authors (e.g.Thuburn and Haine, 2001;Vukicevic et al., 2001) have discussed the undesirable properties of the discrete adjoints to such schemes.A number of studies (Hourdin et al., 2006;Henze et al., 2007;Hakami et al., 2007;Gou and Sandu, 2011) have made calculations showing that the exact adjoint to a forward scheme subject to advective nonlinearities produces unphysical sensitivities and show that more physically reasonable results are obtained by FDA methods such as the Eulerian backtracking method adopted here.One issue is that, when flux limiters are on, the associated nonlinearity is so strong that any tangent linear model diverges rapidly from the full nonlinear model for perturbations of amplitude relevant to applications.The situation is particularly extreme at zero concentration.In this case no unique tangent linear model exists, as direct sensitivity calculations can give different results for positive and negative perturbations of any size (see e.g.Fig. 2b of Thuburn and Haine, 2001, and surrounding discussion).
Based on the above, our view is that it is near-impossible in the presence of advective nonlinearities such as flux limiters, to obtain adjoint sensitivities that are both accurate (in the sense defined in Sect. 1) and practical (in that they do not contain spurious and unphysical sensitivities due to the scheme's nonlinearity).However, for numerous problems in atmospheric chemistry flux limiters are not needed, notably those in which the key species under investigation has a significant background concentration (e.g.O 3 CH 4 , N 2 O, CO, CO 2 etc.), and their localized sources are not too strong.
In process studies with simple chemistry, it may also be preferable to switch off flux limiters to improve accuracy, provided that small regions of negative mixing ratio can be tolerated for short times.This approach also has the advantage of preserving tracer-tracer correlations (Thuburn and McIntyre, 1997).Below, however, we record some tests to determine the magnitude of the inaccuracies introduced by flux limiters.The focus is on the "worst-case scenario" for the spatial structure of the source, when it is isolated to a single grid cell and the background concentration is set to zero.Much better performance can be expected for smoothed sources and non-zero backgrounds.
Figure 4 compares direct sensitivities for test cases I and II for forward model integrations with and without flux limiters (the left-hand plots are identical to Figs. 1b and 2b respectively).The two centre plots give direct sensitivities calculated from nonlinear forward runs subject to flux limiting.The right-hand plots give the discrepancy between these two cases: with the DL subscript indicating a direct sensitivity calculated with flux limiting.For test case I the maximum and average relative errors are 0.13 and 0.035 respectively.For test case II they are 0.20 and 0.037 respectively.For the longer experiment, test case III, the red squares in Fig. 3b show direct sensitivities calculated with flux limiting.The maximum and average errors relative to the direct sensitivities calculated without flux limiting are 0.021 and 0.012 respectively.The maximum error in particular is significantly smaller than for the shorter test cases.
While these errors are significant, it is to be noted that they are the difference between calculating forward sensitivities with and without flux limiting, the latter of which is linear and can be efficiently replicated by an adjoint calculation and the former of which is not.Since our adjoint is equivalent to the forward scheme run without nonlinear flux limiting this is equally well viewed as a comment on the (significant) effect of flux limiting upon direct sensitivity calculations.These errors represent the sensitivity of the CTM to the use of flux limiting and perhaps reflect the artificial nature of their use to control a concentrated source of tracer in a single grid box given the large spatial gradients that result.We also note that the errors due to flux limiting are significantly smaller for the longer run, in contrast to the error due to density inconsistency which remains of a comparable size for all run lengths.Turning on flux limiting in runs of RETRO-TOM results in an average error of 0.007 for test case II, which is significantly smaller than the 0.037 reported for flux limiting applied to the direct sensitivity calculation.We believe that this decreased impact is at least partially explained by the larger source region (a 10 × 10 × 4 receptor) for the backwards calculation.

Resolution study
So far we have presented results obtained at a horizontal resolution of approximately 2.8 • × 2.8 • .In this section we discuss results obtained when RETRO-TOM is run as an adjoint to TOMCAT with an increased 320 × 160 (≈ 1.1 • × 1.1 • ) horizontal resolution.The new grid has an unchanged vertical resolution and 2.5 times as many grid points in both longitude and latitude.As a consequence of the increased horizontal resolution the advective time step is reduced to 5 min, a reduction by a factor of 6 in comparison to the results previously presented.In total therefore high-resolution runs take approximately 38 times longer than their 2.8 • counterparts, a potentially prohibitive increase for long studies with a large number of tracer species.
To examine the robustness of our results to changes in resolution we recalculate C * for test cases I and II with RETRO-TOM run at 1.1 • resolution.Figure 5 contrasts the sensitivities obtained at 2.8 • and 1.1 • resolutions for both test cases I and II.The 2.8 • results are as in Figs.1a and 2a but have been replotted here using a slightly different colour scale (in order to accommodate the 1.1 • results) and without indicating individual grid cells.Clearly the higher-resolution results resolve finer-scale structures; however, the overall picture remains similar.In order to quantify this effect we consider the globally averaged sensitivity for each of the maps.For test case I we find that the average sensitivity of the highresolution results is just over 30 % higher than it is at lower resolution.For test case II the average sensitivity is just under 20 % higher.For both test cases the higher-resolution results were compared against direct sensitivity results over a small (6-grid-cell) forward patch.As for the 2.8 • results reported above, the relative errors are of order 10 −8 .

Conclusions
In this technical note the development of RETRO-TOM, an adjoint to the "dynamical core" of the chemistry transport model TOMCAT has been presented.RETRO-TOM has been shown to combine the conceptual and numerical advantages of a "finite difference of adjoint" (FDA) model with the accuracy of an "adjoint of finite difference" (AFD) model.The three key aspects of the model development are: 1. the use of Hourdin and Talagrand (2006)'s Eulerian backtracking framework; 2. identification of the "time-symmetry" (Hourdin et al., 2006) of the Prather (1986) advection scheme; 3. the recognition of and careful treatment of the "density inconsistency problem" (Jöckel et al., 2001).
When flux limiters in the advection scheme are not used, for example as is appropriate for the inverse modelling of species with a high background concentration (e.g.CO 2 , N 2 O, CO or CH 4 ), accuracies of up to 10 −8 for a 6-month transport problem are recorded.With flux limiters on, there are good arguments, supported by calculations (see Hourdin et al., 2006, and above), for saying that it is better to use the Eulerian backtracking model formulated here than to attempt to create an exact adjoint to the nonlinear scheme.
One advantage of a CTM adjoint formulated in the Eulerian backtracking framework is that it can be used for process studies designed to investigate the origin of air masses, in exact analogy with Lagrangian back trajectory studies (e.g.Manning et al., 2003;Stohl et al., 2003;Seibert and Frank, 2004;Stohl, 2006;Polson et al., 2011).By way of demonstration, RETRO-TOM has recently been used successfully (Haines and Esler, 2014) to quantify the efficiency of different surface source locations for very-short-lived species of halogenated hydrocarbons, in terms of their potential contribution to the halogen flux into the stratosphere.
RETRO-TOM will be suitable, in combination with an adjoint to the tangent linear scheme of TOMCAT, for the various applications (data assimilation, inverse modelling, parameter sensitivity estimation) detailed in the introduction, where it is to be hoped that the higher level of accuracy can be used to exclude numerical errors as a significant factor in some problems, and will prove particularly advantageous in problems where consistency between forward and adjoint models is of vital importance.
A recent development is that an alternative adjoint model for TOMCAT (ATOMCAT, with associated nonlinear inverse model INVICAT), based on the AFD framework, has been under concurrent development by the Leeds University group (Wilson et al., 2013).The ATOMCAT model has an advantage over RETRO-TOM in that it is coupled to the INVICAT model, which is designed to solve nonlinear inverse problems by iterating forward TOMCAT and adjoint ATOMCAT calculations.As a result, at the present stage of development ATOMCAT/INVICAT are suited to a wider range of applications compared to RETRO-TOM.However, ATOMCAT also has certain disadvantages relative to RETRO-TOM including the disadvantages of AFD numerical schemes detailed above, which may be particularly severe when flux limiters are in use (see discussion above).It is also likely to be the case that RETRO-TOM is much better suited as an alternative to Lagrangian backtracking (e.g.Haines and Esler, 2014), because ATOMCAT requires output from a forward calculation at every time step, making it difficult to use in problems that are formulated without reference to a forward calculation.Finally, RETRO-TOM has the advantage that, unlike ATOMCAT which relies on code generated by differentiation of TOMCAT's forward code, RETRO-TOM utilizes the machinery of TOMCAT itself (e.g.Prather advection scheme, column matrix approach to parameterizations), making it relatively straightforward for users to adapt and update in conjunction with the forward model.An interesting topic for future study is the question of whether the ATOMCAT "dynamical core" can be replaced by that of RETRO-TOM without any degradation in model performance, thus creating a full inverse model with the advantages of both schemes.

P. E. Haines et al.: The RETRO-TOM adjoint model
Appendix A: Time symmetry of the Prather scheme A1 First-order moments Following Prather (1986), define orthogonal first-order polynomial basis functions on the interval X 1 < x < X 2 , Note that each function has been normalized so that their extrema are invariant under a rescaling of the interval length.
An appropriate choice of moment coefficients s 0 and s x then allows us to approximate any function φ(x) defined on the interval as φ(x) ≈ s 0 K 0 +s x K x .Here φ is assumed to represent the mixing ratio of tracer within a "grid box" X 1 < x < X 2 of constant air density ρ and so s 0 gives the mean value of the mixing ratio in the grid box.Prather (1986) presented his scheme in terms of moments obtained from integrals of φ over the grid box.For a first-order scheme in one dimension these moments are equivalent to where we explicitly allow for the air density ρ to vary between grid boxes and so have box mass m = ρ(X 2 − X 1 ) in place of Prather's volume of integration V .The "zeroth-order moment" S 0 gives the total mass of tracer in the grid box.
However, here we choose to work entirely in terms of the coefficients (s 0 , s x ) and will use the terms moment and moment coefficient interchangeably.
The process of advection by first-order moments involves first splitting each grid box into several sub-boxes for which appropriate moment coefficients are obtained and then recombining these coefficients in a way that captures transport to neighbouring boxes.Suppose first that we wish to obtain an approximation to φ on a second interval [X 3 , X 4 ] with corresponding basis functions K0 , Kx defined in the same way as Eq. ( A1).An identical polynomial representation φ(x) ≈ s0 K0 + sx Kx is obtained in terms of the new basis functions by taking Now let functions φ (1) , φ (2) and φ (3) denote the mixing ratio in boxes with mass m (1) , m (2) and m (3) respectively, for which as before.We wish to combine these three boxes into a new box of length X weighting the length of each sub-box according to the proportion of the total mass it contributes.Thus the moments from sub-box 1 will take up the portion of the box 0 < x < X1 , that from sub-box 2 X1 < x < X2 and that from sub-box 3 X2 < x < X, where and m = m (1) + m (2) + m (3) is the total mass of the new box.A key feature of the chosen basis functions is that we can change the interval that a function φ (j ) varies over just by changing the basis functions, that is we do not have to recalculate the coefficients s 0 , s x .Thus the moments from each sub box can be combined to obtain a new function, where s (j ) 0 and s (j ) x are as before but K (i) 0 and K (i) x now refer to basis functions defined on the new subintervals.We then construct an approximation φ = s0 K 0 + sx K x to φ on the new interval such that giving s0 = α (1) s (1)

A2 Advection by first-order moments
We now consider one dimensional advection as depicted in figure 6 with moments s 0,i and s x,i defining the mixing ratio profile φi (x) for grid box i.Note that the mean box mixing ratio s 0,i is equal to c i as defined in Section 3.1, and similarly S 0,i , the integral Eq. (A2) involving φi over the grid box, is equivalent to S i .Note also that, whereas the other transport operators that comprise TOMCAT were introduced in Sect.3.1 in terms of their action upon the tracer mass, S, it is more convenient to describe the advective transport in terms of the coefficients s 0,i and s x,i .We start by dividing each grid box into three parts.For box i this gives a left subinterval [0, X L ] representing the possible flux from box i to box i − 1, a right subinterval [X R , X] representing the possible flux from box i to box i + 1 and a central part [X L , X R ] that remains in box i.The lengths of the subintervals are chosen in accordance with the fraction of the total box mass, m i , in each sub-box, where L i is the mass flux to box i − 1 and R i the mass flux to box i + 1.Note that both R i−1 and L i represent a flux between boxes i − 1 and i.In most schemes only a net flux is considered and at least one of these will be zero.We choose to include both here since this allows us to construct the scheme for a general mass flux stencil, treating all grid boxes via a single mechanism whether they are subject to convergent, divergent or uni-directional winds.This has the side effect of generating a scheme general enough to handle bi-directional exchange between grid boxes, a feature which could possibly prove advantageous in modelling a diffusive process.Applying Eq. (A3) we obtain new moments for each sub-box: Once each box has been split in this way we must then recombine the moments appropriately.We need to construct a new box i from the three sub-boxes representing a flux R i−1 from box i −1, a mass m i −L i −R i that has remained in box i and a flux L i+1 from box i + 1. Substituting s (1) 0 = s L 0,i+1 and similarly for s (j ) x into Eq.(A8) and then making use of Eq. (A10) we obtain and and the new box has mass mi Equations (A11) and (A2) allow new moments, (s 0 , sx ), to be obtained from old in a single step.

A3 A matrix representation of advection
For a system of N grid boxes we define a state vector s = (s 0,1 , s x,1 , . . ., s 0,N , s x,N ) T comprised of the moment coefficients associated with each grid box.This is the first-order, one-dimensional, incarnation of the extended state vector s introduced in Sect.3.2.Here it has two components per grid box; for second-order, three-dimensional advection as implemented in TOMCAT each grid box requires 10 components.We can express the process of advection in terms of a 2N × 2N matrix operator A acting on the vector s: with the entries of the matrix A determined from Eqs. (A11) and (A2).We focus on the entries in rows 2i − 1 and 2i of A and write in order that the entries may be given more compactly in terms of 2 × 2 sub-matrices X i , Y i and Z i .Defining T i = s 0,i , s x,i T and T i = s0,i , sx,i T to be the moment coefficients in box i at the start and end of the step respectively we have with X i giving the contribution of the moments in box i − 1 to box i due to a rightwards mass flux, Z i the contribution of the moments in box i + 1 to box i due to a leftwards mass flux and Y i the contribution of tracer that remains in box i.
The entries of X i , Y i and Z i are obtained from Eqs. (A11) and (A2) and are as follows: where the substitutions i − α To apply the Prather scheme to Eulerian backtracking the direction of the mass fluxes must be reversed, together with the box masses at the start and end of each step: .
which can be shown to hold from Eqs. (A17)-( A19) and (A24)-(A26) upon making use of the definitions of the various α and the relation This completes the demonstration of the time symmetry of Prather's first-order moments advection scheme in one dimension.

A5 Extension to second-order moments and higher dimensions
The procedure for second-order moments in one dimension proceeds by direct analogy with that for first-order moments, but with a considerably larger amount of tedious algebra.For second-order moments in one dimension there are three moment coefficients T = (s 0 , s x , s xx ) per grid box, the last of which is associated with the second-order polynomial basis function The matrices X i , Y i and Z i are now of size 3 × 3. We list the nine components of Y i (with the subscript i dropped to avoid clutter) as an illustration: Note that the four entries in the top left corner of Y i are identical to those given previously for first-order moments (A18).However, the five new entries are of greater complexity than their first-order counterparts.We have introduced the following additional definitions to save space: 2α(2) i − α (1) i − 1).(A38) Similar expressions have been obtained for X i and Z i .
The weighing matrices W − and W + for second-order moments in one dimension take the same form as in Eq. (A32) but are of total size 3N × 3N , with the N sub-matrices W i now of the form The key point of our presentation in Eqs.(A1)-( A4) is that the definition of time symmetry, given by Eq. (A33), and its local reduction (A34), remains unchanged for this extended system.To show time symmetry, it is necessary only to calculate X * i , Y * i and Z * i , and then to verify the expressions in Eq. (A34).For example, Y * i can be constructed from Eq. (A37) by making use of the relations (A21) (note that λ * i = λi ).Extension to three dimensions requires storage of the 10 coefficients T = (s 0 , s x , s y , s z , s xx , s yy , s zz , s xy , s xz , s yz ) T (A40) for every grid box.The associated basis functions can be obtained by analogy with Eqs.(A1) and (A36) together with non-symmetric permutations of K xy = K x K y .Higher-dimensional advection is performed in a single dimension at a time, treating the coefficients associated with variation in the other dimensions in terms of their variation in the dimension at hand.For example, for advection in the x direction the coefficients that relate to no variation in y and z (s 0 , s x , s xx ) are treated as for a one-dimensional second-order scheme, those that relate to linear variation in y (s y , s xy ) or in z (s z , s xz ) are treated as for two separate one-dimensional first-order schemes in x and the remaining coefficients s yy , s yz and s zz are individually advected as for three zeroth order schemes.If the moment coefficients are ordered appropriately the associated 10 × 10 matrices X i , Y i and Z i are of block diagonal form with one 3 × 3 block, two 2 × 2 blocks and three diagonal entries.Time symmetry then follows from the time symmetry of each of the sub-blocks.

Figure 1 .
Figure 1.Comparison of sensitivities C * for test case I. (a) Sensitivity calculated from a single run of RETRO-TOM by Eulerian backtracking.The white rectangle indicates the location of the receptor region and the red rectangle illustrates the location of the source patch used in subsequent subplots.(b) The same sensitivity calculated directly from a 10 × 10 patch of TOMCAT forward runs.Results are indistinguishable from those shown in (a).(c) Relative error of the Eulerian backtracking results in (a) in comparison to the (true) direct sensitivity in (b).(d) As for (c) but with Eulerian backtracking performed without reversing the order of the operators.(e) As for (c) but with Eulerian backtracking performed by naive RETRO-TOM, that is without initializing the density field to agree with that in the forward problem.Note the logarithmic colour scale for the error plots.

Figure 2 .
Figure 2. Comparison of sensitivities C * for test case II.(a) Sensitivity calculated from a single run of RETRO-TOM by Eulerian backtracking.The white rectangle indicates the location of the receptor region and the red rectangle illustrates the location of the source patch used in subsequent subplots.(b) The same sensitivity calculated directly from a 10 × 10 patch of TOMCAT forward runs.Results are indistinguishable from those shown in (a).(c) Relative error of the Eulerian backtracking results in (a) in comparison to the (true) direct sensitivity in (b).(d) As for (c) but with Eulerian backtracking performed without reversing the convective mass correction.(e) As for (c) but with Eulerian backtracking performed by naive RETRO-TOM, that is without initializing the density field to agree with that in the forward problem.Note the logarithmic colour scale for the error plots.

Figure 3 .
Figure 3.Comparison of sensitivities for test case III.(a) Sensitivity calculated from a single run of RETRO-TOM by Eulerian backtracking; note the logarithmic colour scale.As before, the white rectangle indicates the location of the receptor region.(b) Sensitivity against longitude along the red line shown in (a) (φ ≈ 38 • S).The solid black line gives C * obtained from RETRO-TOM as previously reported in (a).The points indicate sensitivities obtained from multiple forward model runs performed both with (red squares) and without (blue triangles) flux limiting.

Figure 4 .
Figure 4. Comparison of direct sensitivities calculated with and without nonlinear flux limiting for test cases I (top) and II (bottom).The left hand plots repeat the direct sensitivity results previously shown in Figs.1b and 2b whereas the centre plots show the same direct sensitivity calculated with flux limiters turned on.Note that the colour scales for the four left-hand plots are as used in Figs.1ab (top) and 2a-b (bottom).The colour scales shown here are for the right-hand plots, which give the size of the relative error between the direct sensitivity calculated with and without flux limiting.

Figure 6 .
Figure 6.Schematic showing the incoming and outgoing mass fluxes for box i together with the proportions of the box that they occupy.
. The state vector s * containing moment coefficients associated with Eulerian backtracking is calculated from s * according tos * = Bs * (A22)with the entries of B obtained from those of A by replacing all variables by their Eulerian backtracking (starred) equivalents.Here we focus on the non-zero entries of B in columns 2i − 1 and 2i: indicating that all variables in the 2 × 2 submatrices Eqs.(A17)-(A19) are to be replaced by their Eulerian backtracking equivalents in accordance with Eq. (A21):