MALTA: A Zonally Averaged Global Atmospheric Transport Model for Long‐Lived Trace Gases

We present a two‐dimensional, zonally averaged global model of atmospheric transport named MALTA: Model of Averaged in Longitude Transport in the Atmosphere. It aims to be accessible to a broad community of users, with the primary function of quantifying emissions of greenhouse gases and ozone depleting substances. The model transport is derived from meteorological reanalysis data and flux‐gradient experiments using a three‐dimensional transport model. Atmospheric sinks are prescribed loss frequency fields. The zonally averaged model simulates important large‐scale transport features such as the influence on trace gas concentrations of the quasi‐biennial oscillation and variations in inter‐hemispheric transport rates. Stratosphere‐troposphere exchange is comparable to a three‐dimensional model and inter‐hemispheric transport is faster by up to 0.3 years than typical transport times of three‐dimensional models, depending on the metric used. Validation of the model shows that it can estimate emissions of CFC‐11 from an incorrect a priori emissions field well using three‐dimensional (3D) mole fraction fields generated using a different 3D model than which the flux gradient relationships were derived. The model is open source and is expected to be applicable to a wide range of studies requiring a fast, simple model of atmospheric transport and chemical processes for estimating associated emissions or mole fractions.


Introduction
Supplying observation-based global emission estimates for long-lived gases is important for assessing the success of international policies designed to minimize these emissions and their environmental impacts.A prime example is the quadrennial Scientific Assessment of Ozone Depletion (World Meteorological Organization (WMO), 2022), which reports on emissions of ozone-depleting substances and other greenhouse gases of interest to the Montreal Protocol on Substances that Deplete the Ozone Layer.The Montreal Protocol controls production and consumption of ozone-depleting substances, and other greenhouse gases through the Kigali Amendment, and has been the driver behind ozone layer recovery and the avoidance of large-scale additional global warming (Newman et al., 2009).Global emissions for the majority of greenhouse gases and ozone depleting substances reported in the last Scientific Assessment of Ozone Depletion were quantified using a well-established method, which uses measurements of atmospheric concentration near the Earth's surface coupled with a simple box model of global atmospheric transport in an inverse framework (e.g., Cunnold et al., 1983;Rigby et al., 2013).Here, we limit the term box model to simple low-resolution (1-12 regions for the entire globe) paramaterisations of zonally averaged transport.Simple models of atmospheric transport and loss, such as in a box model representation of the global atmosphere with annually repeating dynamics, have limitations.For example, features of large-scale dynamics, such as the quasi-biennial oscillation of the direction of the zonal wind in the tropical lower stratosphere (QBO), affect surface concentrations of trace gases, but cannot be represented in those models (Montzka et al., 2021;Ray et al., 2020;Ruiz et al., 2021).Low spatial resolution in simple box models also leads to substantial representation errors, when point measurements are compared to simulations of concentrations in very large boxes.Such errors in modeled surface concentration are propagated into emissions estimates if the measurements used are not sufficiently representative of the spatial model extent that they represent, when those emissions are derived from measured mole fractions.
There is a need to improve the interpretation of temporal changes in atmospheric concentration measurements, to supply more accurate estimates of emissions on policy-relevant timescales, such as annually.For example, emissions of CFC-11, quantified using a box model with annually repeating mixing and advection, increased after their global phase-out (Montzka et al., 2018), most likely due to production that was not reported to the Montreal Protocol (Park et al., 2021).The initial box model emissions estimates were adjusted to account for variability in atmospheric dynamics using three-dimensional (3D) models, revealing large uncertainties in the derived emissions for some years when annually repeating transport was used (Montzka et al., 2021).However, this procedure of using a 3D model to correct for variability in atmospheric dynamics is complex and not readily available for other species due to the computational and personnel expense of the procedure.Therefore, to account for interannual variability in atmospheric dynamics, and reduce representation errors, there is a need to improve the modeling of atmospheric trace gases.Simultaneously, any transport modeling approach must be fast and simple to use, given the number of gases that must be quantified for reports such as the Scientific Assessment of Ozone Depletion, from multiple measurement networks, and interim scientific studies (e.g., Mühle et al., 2022;Stanley et al., 2020;Western et al., 2023).This makes quantifying emissions using established three-dimensional chemical transport models prohibitive, except for individual focused studies (e.g., Claxton et al., 2020;Lunt et al., 2015;Montzka et al., 2021).
Zonally averaged two-dimensional (2D) models of atmospheric transport have been widely used to study problems in the atmospheric sciences (Fleming et al., 2020;Hough, 1989;Strand & Hov, 1993).Many past studies using 2D models have focused on processes in the middle atmosphere (∼10-90 km).A 2D model with annually repeating transport, which is not openly available to users, has previously been used to quantify emissions of long-lived halocarbons (e.g., Adcock et al., 2018;Laube et al., 2016;Newland et al., 2013).
To address the need for a medium complexity modeling approach, here we introduce a new zonally averaged twodimensional model of global atmospheric transport, named MALTA: Model of Averaged in Longitude Transport in the Atmosphere.MALTA uses more representative transport, at higher resolution in space and time, than current box-modeling approaches, and considers time-varying large-scale dynamics.Unlike previous 2D models, MALTA is open-source and aims to be accessible to a broad community of users, with the primary function of quantifying emissions of greenhouse gases and ozone depleting substances with improved representation over current box models.The model and its components are written in the widely used Python programming language, which should help accessibility, compared to compiled languages.Performance is greatly improved over a standard Python implementation by exploiting just-in-time compilation (Lam et al., 2015).Section 2 details the derivation of the model transport parameters; Section 3 describes the numerical schemes applied to solve model transport, sources and sinks; Section 4 presents some approaches to assess the performance of the model; and the concluding remarks are presented in Section 5.

Derivation of Zonally Averaged Transport Parameters
Zonally averaged transport characteristics of the atmosphere are paramaterised from three-dimensional meteorological fields and transport tendencies (see e.g., Plumb & Mahlman, 1987) for each calendar month from 1980 onward.Prior to 1980, MALTA uses a monthly climatology  of the transport parameters.The model is in Cartesian coordinates, where y is the horizontal meridional coordinate of meters from the equator, defined as where R is the radius of Earth and ϕ is latitude in radians.The vertical component z is a scale-height where p is zonally averaged atmospheric pressure, H = 7,200 km and p s = 1,000 hPa.The horizontal (y) resolution is every 10°in latitude and the bounds of altitude (z) are at 30 levels (29 layers) spaced equally in logpressure coordinates between 1,000 and 10 hPa (approximately 1,150 m).The resolution was chosen to attempt a balance between representation and fast computation.Although not used here, unequal vertical spacing would also be valid.
No topography is present in MALTA, and thus there is no explicit parameterization of phenomena such as topography-induced gravity waves.The choice of Cartesian coordinates has been made for their relative simplicity over more complex, but perhaps better surface-resolving, coordinate systems (see, e.g., Tung, 1982;Iwasaki, 1989).
By denoting a zonal-temporal average by an overbar and a prime as the deviation from this average, the zonally averaged budget equation for a trace gas with an atmospheric mole fraction q is, where v is the north-south wind velocity and w is the upward wind velocity (we define u = (v,w)) , and s represents the sources and sinks.The derivative ∇ = ∂ ∂y ĵ + ∂ ∂z k, where ĵ and k are unit vectors.Here u′q′ is an eddyflux term, which we expand on in the following section.

The Zonally Averaged Flux-Gradient Relationship
The eddy-flux term u′q′ from Equation 3 can be paramterised as where K is a 2 × 2 eddy transport tensor, The eddy transport tensor can be expressed as the sum of a symmetric diffusive tensor, D, and antisymmetric advective tensor, A, (Bachman et al., 2020;Griffies, 1998;Plumb, 1979), This allows the flux-gradient relationship in Equation 4 to be rewritten, where ψ is the streamfunction for zonally averaged coordinates and Following this approach, the eddy-induced advection v* and w* is derived (Griffies, 1998) following and the eddy-induced diffusion terms are taken directly from the tensor components of D from Equation 5, that is, Substituting this flux-gradient relationship into Equation 3 gives the continuity equation in terms of a residual velocity and eddy-induced diffusion, Alternate derivations of residual circulation exist (see e.g., Andrews et al., 1999;Holton, 1981).

Tracer-Based Estimates of Eddy Transport
To estimate the components of K we follow the approach of Bachman et al. (2015).The idea is to use the eddyflux term in Equation 4 to estimate Tracer experiments from a three-dimensional model provide the information to diagnose 2D transport (Plumb & Mahlman, 1987).This approach has been used, at least in part, to derive transport parameterizations in various 2D models (e.g., Fleming et al., 2020;Hough, 1989;Strand & Hov, 1993).Alternatively, residual circulation has been derived directly from reanalysis fields, separately from the tracer-derived diffusion fields (see e.g., Andrews et al., 1999;Holton, 1981, for this derivation.).
We choose to use six sinusoidally orthogonal pairs of tracer fields (varying in y and log(z)) to initialize the model (following Bachman & Fox-Kemper, 2013, Appendix B).The tracers are inert (with a molar mass nominally equal to carbon monoxide) and have no emission source.The experiments are run using the GEOS-Chem transport model (Bey et al., 2001; The International GEOS-Chem User Community, 2021) driven by MERRA-2 reanalysis meteorology (Gelaro et al., 2017) at 4°× 5°spatial resolution from 1980 onward.Improved spatial resolution of 3D models is known to generally improve representation of trace gas distributions (e.g, Stanevich et al., 2020).Improved transport in the 3D model would yield improved derived zonally averaged transport, however we deem 4°× 5°spatial resolution sufficient given the coarser resolution of MALTA.Each tracer field is run for 1 month from the initialized state to infer the mean zonally averaged eddy-transport for that month, and the three-dimensional tracer and meteorological fields are output every hour.The three-dimensional fields are used to derive the eddy transport tensor.
We derive the eddy-transport by first interpolating the three-dimensional fields onto the two-dimensional vertical and latitudinal resolution at each time step using bilinear interpolation (Zhuang et al., 2020).Gradients are calculated using second-order central differences.The choice of interpolation scheme will introduce small differences to the calculated gradients but have little impact on the resultant transport.We infer the tensor K using the tensor of tracer fields q using an approach akin to maximum likelihood estimation, where i and j are an index of the vertical and horizontal grid cells, by and The square matrix P ij is a diagonal weighting matrix, where 1 is a column vector of ones (similar to, e.g., Bachman et al., 2015;Bratseth, 1998).
Some quality control, following Plumb and Mahlman (1987), is placed on the derived diffusion parameters.That is, D yy (which is equal to K yy ) was set to a minimum value of 10 4 cos 2 ϕ m 2 s 1 , and we ensure that D 2 yz ≤ D zz D yy (ensure antidiffusion is not possible) by adjusting D zz accordingly.

Boundary Layer Transport
Mixing timescales in the lowest model levels (the boundary layer) are rapid and, in a zonally averaged representation, transport may intersect with topography.Problems with using the outlined approach in Section 2.2 to derive flux-gradient relationships in the lower model levels have been previously identified (Plumb & Mahlman, 1987).This is part of a wider problem with representing a zonal average near the surface in a consistent coordinate system.In this work, we have made adjustments to the transport in the lowest two layers of the model to compensate for these shortcomings.We retain the advective fields as per Section 2.3 and meridonial diffusion, and adjust the vertical diffusion to better represent the zonal transport through optimisation.First, we assume that diffusion is controlled only by the vertical and horizontal components (D yz = 0).The vertical diffusion coefficients are subsequently scaled to best match that of the surface concentrations from a 3D model run, using the set up described in Section 4.1 to produce surface mole fractions, and the goodness of fit is evaluated using the root-mean-square error.The diffusion coefficients are each scaled from their initial values derived in Section 2.2, and the scaling is applied to each year and month equally.The scaling is applied as a latitudinal gradient, where the scaling coefficients are optimized at the poles and the equator and the scaling is linearly interpolated between these values.Optimisation is performed using brute-force optimisation.

Advective Fields
The advective fields are taken as zonal and monthly averages from MERRA-2 reanalysis for each year, and combined with the eddy-induced residual advection (Section 2.1) to create a single advective field.Following interpolation and averaging, the 2D advective fields become divergent.We address this problem by adjusting the vertical velocity w + w * following Shine (1989) to ensure that the area-averaged vertical velocity is zero.That is, in each layer i, w i + w * i is additively adjusted by a factor of Using the adjusted values for w + w * , we compute the non-divergent values for v + v * by ensuring mass balance during advection between grid cells.The adjustment to w + w * is typically ∼2% and to v + v * is typically ∼10%.

Computation of Model Processes
The parameters derived in Section 2 transport the mixing ratio of the trace gas following Equation 9.In addition to the transport in Section 2, a convective parameterization has been implemented to approximate the zonally averaged convection.The transport, source and sink processes are broken down using time splitting.Initially, at time t, the tracer is emitted into the atmosphere (generally in the lowest layer of the model).Following this, at time t + τ 1 , the field is advected using the residual velocity (Section 3.1).At time t + τ 2 , diffusion occurs in the y and z direction using the diagonal of the eddy-diffusion tensor (using D yy and D zz , Section 3.2), followed by the offdiagonal eddy diffusion component at time t + τ 3 (Section 3.2.1).Convection occurs at time t + τ 4 (Section 3.3).Finally at time t + τ 5 , any atmospheric sink is imposed (Section 3.4).The numerical implementation of these approaches is described in the following subsections.

Advection
The advection scheme is based on Lin and Rood (1996).A brief overview is given here; see the Supporting Information S1 for a more detailed description of the steps and equations used.The scheme conserves mass, a uniform concentration field, is positive definite, ensures monoticity, is second-order accurate and is numerically efficient.The tracer field is advected using the residual velocities v + v * and w + w * and computed on an Arakawa C-grid (Arakawa & Lamb, 1977).A monoticity constraint on the upwind fluxes is imposed following Equation 5in Lin et al. (1994), meaning that a monotonic mixing ratio gradient before advection will maintain this monotonic distribution after advection.The fractional fluxes in the Lin and Rood (1996) scheme are computed using the Piecewise Parabolic Method (PPM), (Colella & Woodward, 1984).Note that, given the parameters in Section 2 and a time step of around 8 hr, the Courant number does not exceed one and so the semi-Lagrangian component presented in Lin and Rood (1996) is not needed.Over and undershoots are eliminated using the first constraint proposed in Appendix C of Lin and Rood (1996), where edge values are computed following Carpenter et al. (1990).Finally, fluxes between grid cells are computed following Colella and Woodward (1984).

Diffusion
Diffusion is controlled through the diffusion tensor D, Section 2.1.We solve the diffusion due to the diagonal of the eddy-diffusion tensor (D yy and D zz ) using an Eulerian approach on an Arakawa C-grid (Arakawa & Lamb, 1977).Diffusion is calculated using second-order central finite difference in the vertical and horizontal.We use a fourth order Runge-Kutta, or RK4, method for the temporal discretization (e.g., Butcher, 1996).
While finite difference methods are sufficient for diffusion from the diagonal of the eddy diffusion tensor, the offdiagonal elements need more detailed treatment.

Mixed Derivative Diffusion
Mixed derivative second order finite difference methods are not generally positivity-preserving, which is necessary for the transport of atmospheric trace gases.Instead, we adapt the approach proposed by du Toit et al. ( 2018) to ensure positive definite tracer concentrations following the diffusion terms containing a mixed derivative.The diagonal diffusion terms (using D zz and D yy ) are always computed using the Eulerian scheme first described in this section.
The idea is to reframe the mixed-derivative diffusion processes as advective processes that can be solved using a positivity-preserving scheme, as in Section 3.1.Diffusion using the D yz component can be written where and Equation 15 takes a similar form to the advective terms in Equation 9.The ∂q ∂z term in Equation 16 can be solved using standard second order central finite difference methods.However, following Equation 16, as q → 0,v d → ∞.We therefore impose the limit that v d = 0 when q < ϵ, where ϵ is around five orders of magnitude less than the expected surface concentration of the tracer.Conceptually, this implies that the diffusion tensor is diagonal for small concentrations.The same approach is taken for the corresponding w d for diffusion using the D zy component.
The terms v d and w d are calculated using an Arakawa A-grid (Arakawa & Lamb, 1977), but the advection scheme in Section 3.1 requires an Arakawa C-grid.The transport scheme therefore requires interpolation of v d and w d onto the Arakawa C-grid using linear interpolation.These advective-like terms are then used as described in Section 3.1 to advect, and therefore diffuse, the tracer field.

Convection
Zonally averaged transport thus far is generally valid when the timescales of transport is much shorter than the atmospheric lifetime of the atmospheric tracer.The transport of atmospheric tracers with shorter lifetimes may not be adequately simulated due to rapid vertical convective transport, especially in mid-latitudes and the tropics (Chatfield & Crutzen, 1984).We follow the approach of Strand and Hov (1993) to simulate vertical convective transport.Convection is computed as a one dimension column model for each latitude band (similar to, e.g., Pleim & Chang, 1992).
Air is drawn out of the boundary layer (here, we define the boundary layer as ≲2,000 m) and detrained into the atmosphere above the boundary layer.The mixing inside the boundary layer occurs both upward and downward, and subsidence occurs to the level directly below each model level above the boundary layer.The detraining mass flux is the zonal average of the fields taken from MERRA-2 reanalysis, and averaged for each month.Each time step, for each one dimensional column, is solved using a Crank-Nicolson method following Strand and Hov (1993), Appendix E.

Sinks
An explicit chemistry scheme is not currently implemented in the model.Loss due to reaction with the hydroxyl radical (OH) is calculated using the Arrhenius equation, with reaction rates taken from Burkholder et al. (2020) and temperatures taken from the MERRA-2 reanalysis.The OH fields were taken from Pimlott et al. (2022) using their values derived for 2010, with different OH fields (concentration and distribution) applied for each month, which were zonally averaged and scaled, by the same constant value in space and for each month, to the tropospheric lifetime, or lifetime due to OH reactive loss, of methyl chloroform (6.1 years, Burkholder & Hodnebrog, 2022).
Stratospheric sinks due to photolysis and reaction with excited oxygen (O 1 D) are calculated using a first order loss frequency field (the inverse of the local lifetime in each grid cell), taken from zonally averaged loss fields in the climatological stratosphere (driven by the MERRA-2 reanalysis) from prior 3D model simulations (e.g., Murray et al., 2012).We assume that only loss due to OH occurs in the troposphere.Given the tropospheric lifetime, we tune the stratospheric lifetime such that the total atmospheric lifetime is equal to that reported in Burkholder and Hodnebrog (2022), unless another desired lifetime is needed, by summing the inverse lifetimes in the troposphere and stratosphere.The atmospheric lifetime is estimated when the target species is in near-steady state in the atmosphere, following a 20-year spin-up period.These loss frequencies vary by month but repeat annually.Using a tuned loss field may lead to a correctly derived atmospheric lifetime of a species despite errors in transport.However, we justify the use of the offline loss field based on the improved computational speed and ease of implementation.

Model Performance
In this section we evaluate the performance of the model.In Section 4.1 we compare mole fractions of the atmospheric tracer sulfur hexafluoride (SF 6 ) modeled using MALTA to those from GEOS-Chem, which was used to derive the eddy-transport tensors, and a 12-box model.Where we compare atmospheric concentrations at a particular altitude, the scale-height in MALTA is treated as meters above sea-level (masl).Section 4.2 evaluates the timescales of inter-hemispheric transport of MALTA compared to GEOS-Chem and those from other models and observations presented in previous literature, and Section 4.3 evaluates the stratosphere-troposphere exchange, using SF 6 as a tracer.Section 4.4 assesses the ability of MALTA to represent the influence of the QBO on atmospheric mole fractions of CFC-11, and discusses some the observed features.We derive emissions of CFC-11 in Section 4.5 using modeled mole fractions from a distinct 3D transport model with a known emissions input.

Modeled Surface Mole Fractions
Here we compare modeled surface mole fractions of SF 6 taken from MALTA to those output from the 3D model GEOS-Chem.In addition, we compare these surface mole fractions to those modeled using a 12-box model (Cunnold et al., 1983;Rigby et al., 2013).Note that MALTA has been tuned to the transport in GEOS-Chem, whereas the 12-box model has not.The 12-box model mole fractions will differ from MALTA or GEOS-Chem but the inclusion of a simple independent model for completeness is nonetheless informative.We model SF 6 as a tracer with an infinite lifetime (in MALTA, GEOS-Chem and the 12-box model) to assess only the transport within the model.Emissions are from the EDGAR v4.2 emissions inventory (Janssens-Maenhout et al., 2011).Emissions are assumed to be constant for each month in a given calendar year.that the initial conditions in 1980 for the modeled simulations are somewhat higher than those from atmospheric measurements (Simmonds et al., 2020).
The MALTA and GEOS-Chem simulated mole fractions generally agree well when compared to the grid cells in GEOS-Chem containing 12 measurement sites in the National Oceanic and Atmospheric Administration's (NOAA) measurement network, which are detailed in Table S1 in Supporting Information S1. Figure 1 shows a comparison of the modeled monthly mean SF 6 surface mole fraction simulated using MALTA, GEOS-Chem and the 12-box model between 1980 and 2008 at four site locations: NWR, MLO, SMO, and SPO, and the respective growth rate in the mole fractions.When compared to GEOS-Chem, the root-mean-squared error of the monthly mole fraction over all sites for MALTA is 0.05 ppt, compared to 0.11 ppt for the box model.The bias in the monthly mole fraction in MALTA is 0.02 ppt for all sites compared to 0.08 ppt in the 12-box model.MALTA does not outperform the 12-box model at every site, with MALTA slightly overestimating the mole fractions in the Southern Hemisphere compared to GEOS-Chem and the 12-box model, most likely due to the faster interhemispheric transport time in MALTA (see Section 4.2).The average seasonal cycle at each site (see Figure S1 in Supporting Information S1) is generally similar between MALTA and GEOS-Chem, with the exception of SMO.At SMO, the agreement between the 12-box model and GEOS-Chem is better than that between MALTA and GEOS-Chem.The SMO site (in American Samoa) is known to have strong interannual modulation, which can be somewhat different to its zonal average, due to changes in the South Pacific Convergence Zone and other large-scale features (Hartley & Black, 1995;Raiter et al., 2020).
Figure 2 shows the monthly mean mole fraction for January 2005 of the 12-box model, MALTA and the zonal average of GEOS-Chem throughout the atmosphere, demonstrating the difference in their horizontal and vertical resolution and simulated mole fractions.

Inter-Hemispheric Transport
We evaluate the inter-hemispheric transport of GEOS-Chem and MALTA using two metrics: the age of air of SF 6 between the northern and southern hemisphere mean surface mole fractions (Waugh et al., 2013) and the interhemispheric exchange time, τ ex , (Patra et al., 2009) using the model runs described in Section 4.1.We follow a similar approach to Yang et al. (2019) when evaluating these metrics in order to compare to their results.The age of air of SF 6 (a SF 6 , the time lag for the mean mole fraction in the southern hemisphere to equal that in the northern hemispheric) is calculated as where t is time, q is the average mole fraction in the northern and southern hemisphere, denoted by its subscripts.The inter-hemispheric exchange time is Journal of Advances in Modeling Earth Systems where E is the emissions in each hemisphere.As in Yang et al. (2019), the age of air is calculated using 23-month smoothed mole fractions (using a 23-month exponentially weighted moving window) for each site.We use simulated mole fractions from 1995 to 2009, using the grid cells containing the sites BRW, MLO, CGO and SPO (see Table S1 in Supporting Information S1) in the calculation of τ ex ; and MHD, NWR and THD, along with the average mole fraction over 60-90°S, in the calculation of a SF 6 , as Yang et al. (2019).
Note that the modeled mole fractions are not meant to be representative of the truth in this scenario.Instead, we assume that the hemispheric distribution of emissions from EDGAR v4.2 are sufficiently representative of the true emissions distribution to make a meaningful qualitative comparison between MALTA, GEOS-Chem and previously derived metrics from other models and measurements.2019) is slightly shorter (around 1.4 years).The shorter timescale for a SF 6 in MALTA indicates that MATLA may have some shortcomings in simulating transport from northern hemisphere midlatitudes to the northern hemisphere tropics, that is, transport is too rapid (Yang et al., 2019).
The mean annual cycle of τ ex for each month is shown in Figure 3a.The age of air of MALTA is generally younger (exchange is faster) than GEOS-Chem outside of January to May.The seasonal cycle between MALTA and GEOS-Chem share some similar features (e.g., the most rapid transport occurring during the late northern hemisphere summer).However, MATLA does not exhibit the large slowing of interhemispheric transport in the northern hemisphere winter as seen in GEOS-Chem.The reason for the differences between the seasonal cycles in GEOS-Chem and MALTA may come from the zonal averaging of the parameters within the model, for example, where the distribution of the zonal transport is strongly non-normal.The transport processes controlling inter-  hemispheric transport is known to be largely due convective processes, in particularly in tropical and sub-tropical latitudes (Gilliland & Hartley, 1998;Lintner et al., 2004;Patra et al., 2009).The cross-equatorial transport is dominated by large-scale convective outflow (Hartley & Black, 1995).Regional, longitudinally dependent convective processes can modulate interhemispheric transport (Lintner et al., 2004), which are clearly neglected in MALTA given the lack of any longitudinal variations.For example, the Indian Ocean Monsoon, and strengthening of convection along the Indian Ocean Intertropical Convergence Zone, is known to modulate interhemispheric transport during the Northern Hemisphere summer (Lintner et al., 2004).This may indicate a shortcoming in MALTA in representing interhemispheric transport that is dominated by convective process, compared to the advective dominated processes in extra-tropical latitudes (Patra et al., 2009).
Shortcomings in representative convective processes in MALTA may have further implications for modeling of gases that have strong zonal variability and can be rapidly transported into the upper troposphere and lower stratosphere by convective processes.This will be particularly pertinent for short-lived ozone-depleting substances (lifetime <6 months) where large emissions sources are co-located within the Asian monsoon region (e.g., Adcock et al., 2021), and may mean that this model will introduce additional errors when studying these gases.

Stratosphere-Troposphere Exchange
We calculate the stratosphere-troposphere exchange (STE) of SF 6 from 1980 to 2008 using GEOS-Chem and MALTA using the emissions scenario described in Section 4.1.For both models the location of the tropopause is taken from MERRA-2 reanalysis, based on a minimum in the lapse rate, and we define the troposphere as all grid cells below the grid cell containing the tropopause.As SF 6 is modeled as having an infinite lifetime, we define stratosphere-troposphere exchange, F S→T , as where the subscripts S and T denote the stratosphere and troposphere, B is the burden and E is the emissions to the troposphere per unit time.
Figure 4 shows a comparison of the STE of SF 6 from 1980 to 2008 between the two models.Both models exhibit an average STE maximum in June and a minimum in October, as expected (Appenzeller et al., 1996).The mean absolute magnitude of the STE is slightly larger in MALTA than GEOS-Chem (4.2 vs. 3.5 Gg yr 1 ), and GEOS-Chem has, on average, slightly higher transport into the troposphere than MALTA (overall mean STE of 0.7 Gg yr 1 for MALTA vs. 0.8 Gg yr 1 for GEOS-Chem).These differences are very small compared to the amplitude of the seasonal cycle of STE (∼10 Gg yr 1 ).It should be noted that GEOS-Chem is known to underestimate STE compared to some other transport models (Hu et al., 2017).Stratospheric mole fractions of SF 6 in MALTA generally appear to be lower (older) than in GEOS-Chem, despite good agreement in STE between MALTA and GEOS-Chem (Figure 2).The lower stratospheric mole fractions could suggest that stratospheric mixing timescales in MALTA are slower than GEOS-Chem, or perhaps weakened tropical vertical transport in MALTA resulting in older air becoming circulated through the Brewer-Dobson circulation.If tropical vertical transport is too weak, STE may be being compensated through stronger extratropical STE than in GEOS-Chem, giving an overall similar STE.

The Quasi-Biennial Oscillation
A large motivation for the development and use of a 2D model over existing box models is to include features of large-scale dynamics and their interannual variability.The QBO has been shown to influence concentrations of trace gases at the surface, in particular of the long lived trace gas CFC-11 (Montzka et al., 2021;Ray et al., 2020;Ruiz et al., 2021).
We assess the representation of the QBO in MALTA using transport of CFC-11, using an average lifetime of 52 years, tuned from an initial 3D field from Murray et al. (2012).Emissions in MALTA are held constant at 88 Gg yr 1 following a 50 years spin up period (using approximate emissions from Laube et al. (2023)), which provides a near-steady state.
Figure 5 shows the latitudinally averaged partial pressure anomaly of CFC-11 in MALTA as a function of time between 1990 and 2020 (see e.g., Ray et al., 2020).The anomaly has been filtered using a Butterworth bandpass filter of 22-34 months, where the mean QBO period is around 28 months (Baldwin et al., 2001).The QBO-induced partial pressure anomaly can generally be seen to propagate from the stratosphere down to the surface.There is a clear disruption in the signal around 2016, following an anomalous QBO phase (Newman et al., 2016).To a lesser degree, there is some non-continuous propagation of the signal downwards from the stratosphere during the mid-1990s to 2000. Figure 5 shows the east-west equatorial wind anomaly at 20 hPa (filtered in the same way as the partial pressure anomaly), where the winds have been taken from radiosonde soundings from Singapore in order to compare against an independent data set (available as a Supporting data set, Western (2024a)).The phase of the east-west stratospheric wind anomaly at 20 hPa agrees qualitatively well with the phase of the partial pressure anomaly at 20 hPa, which demonstrates that the model appears to qualitatively represent a QBO-induced anomaly well.Analysis of the time-lag between the modeled and observed 20 hPa wind anomaly, using cross correlation, shows a time lag of 2.3 years, or around 28 months, the mean period of the QBO.At 70 hPa, the time lag is 2.1 years (around 25 months), which is within the variability of the period of the QBO.

Estimation of Modeled Emissions
One of the primary foreseen uses of MALTA is to derive emissions of trace gases.Here, we demonstrate this use case by estimating emissions of CFC-11.As emissions are not known exactly in reality, we instead choose to use 3D-model-generated mole fractions, based on an emission history derived from a simple analysis of atmospheric measurements at remote sites across the globe (Montzka et al., 2021).We estimate emissions using modeled mole fraction output from the TOMCAT offline 3D chemical transport model (Chipperfield et al., 2015(Chipperfield et al., , 2016)).The TOMCAT simulation used ERA-5 reanalyses from the European Centre for Medium-Range Weather Forecasts (ECMWF) (Hersbach et al., 2020) and was run at 2.8°× 2.8°resolution at 60 model levels up to 0.1 hPa.The lifetime of CFC-11 in the TOMCAT simulation (54 years) differs from the steadystate lifetime in MALTA (52 years, Section 3.4).See Montzka et al. (2021) for a more detailed description of the TOMCAT set up.The TOMCAT model is not fully independent from GEOS-Chem, as many chemical transport models use similar underlying physics and their parameterization, and both MERRA-2 and ERA5 reanalysis originate from similar assimilation schemes and observations.Nevertheless, we deem TOMCAT a suitable test case simulator.
Emissions used to generate the mole fraction are the smoothed emissions described in Montzka et al. (2021), and are distributed uniformly over the globe.The mole fractions used to estimate emissions were the monthly mean model 3D-output at the location of 12 measurement sites in NOAA's global monitoring network, which spans from 90.0°S to 82.5°N (see Section 4.1; Table S1 in Supporting Information S1).The initial conditions for the mole fraction in MALTA were estimated following a 50 years model spin-up which was then scaled to approximate that of the simulated measurement in the first month minus half the difference in the mole fraction between the first and second month of the TOMCAT output.Emissions were estimated between 2000 and 2019.
Sensitivities in MALTA were generated by perturbing the emissions from an initial reference run by 1 Gg for each latitudinal band for each year.The reference mole fraction was subtracted from the mole fraction output at each measurement site (and divided by 1 Gg) to give the sensitivity of each measurement to emissions.Each forward model run was run from 2000 through 2019.This computation is an embarrassingly parallel problem, where each perturbed forward run can be run independently.This takes approximately 45 min on a four core laptop, and scales for high performance computing depending on the number of available threads.
To estimate annual emissions of CFC-11 from simulated mole fractions, we follow a Bayesian inverse framework, assuming a normal likelihood and prior distribution.Under the assumption of normality, an analytical solution can be found (see e.g., Tarantola, 2005).Under the assumption of normality, we can derive the most probable estimate of the deviation of emissions from the reference run, and the uncertainty in this estimate is Observations are contained in y and sensitivities in H, such that emissions, x, form a linear operation y = Hx + ϵ, where ϵ is the zero mean normal error.Emissions are assumed to be constant in each calendar year (they vary monthly in the TOMCAT simulation) and we infer emissions for each 10°latitude band, which is contained in the vector x.The prior mean emissions are random truncated normal perturbations around the TOMCAT annual emissions, using a standard deviation of 10 Gg yr 1 , where the truncation ensures that only positive prior mean emissions are allowed, and form the vector x a .We then assign a one standard deviation uncertainty of 20 Gg yr 1 on the prior distribution in each latitude band, and form the diagonal of Σ a .The likelihood assumes a combined model-measurement error of 5 ppt on the diagonal of Σ ϵ .Both the likelihood and prior distribution are assumed to be identical and independently distributed.
Figure 6 shows the estimated CFC-11 emissions derived from the TOMCAT simulated mole fractions, the emissions used by TOMCAT to simulate the mole fractions and the prior distribution placed on the emissions.The emissions derived using MALTA show a marked improvement over the prior mean mole fractions.However, greater interannual variability exists in the derived emissions compared to the emissions used to simulate mole fractions in TOMCAT.Likely reasons for the errors are the error in representing 3D fields in two-dimensions, differences in the driving transport and sinks and the lack of spatial information content available in the simulated measurements.

Conclusions
We have presented a 2D zonally averaged global model of atmospheric transport named MALTA, Model of Averaged in Longitude Transport in the Atmosphere.The intention for the model is primarily to simulate longlived greenhouse gases and ozone-depleting substances.As well as halogenated substances, as has been mainly presented, MALTA could be used for other greenhouse gases, such as methane or nitrous oxide, and emissions and transport of non-reactive substances in the atmosphere.A major anticipated application is to infer emissions of these long-lived gases from measurements of atmospheric mole fractions.
The model is driven by monthly averaged transport, derived from MERRA-2 meteorology, and eddy transport derived from experiments using the GEOS-Chem 3D model.The model's zonal transport is generally comparable to GEOS-Chem, although its inter-hemispheric transport appears to be slightly faster than both GEOS-Chem and that derived from measurements.Emissions derived with MALTA from mole fractions simulated using a distinct 3D transport model, TOMCAT, are generally in good agreement with those used as input to the TOMCAT forward simulation.A motivating factor for creating the model was to better simulate the influence of large-scale dynamics on surface mole fractions over existing simplified models, such as box models.MALTA produces a change to surface concentrations and throughout the middle atmosphere, which is in phase with the QBO at 20 hPa.
The model is open-source and freely available.Transport parameters can be derived from different meteorological reanalysis and/or 3D models using an approach similar to that presented.The major intended use is in the estimation of long-lived trace gas emissions, although the investigation of many other processes is possible using the model, such as the study of dynamical transport.

)
Journal of Advances in Modeling Earth Systems 10.1029/2023MS003909 WESTERN ET AL.

Figure 1 .
Figure 1.A comparison of the modeled mole fraction and its growth rate at the grid cells containing four measurement sites using GEOS-Chem, MALTA and a 12-box model.The four represented measurement sites are given by their three letter acronym, and are (a) Niwot Ridge, Colorado, USA, (b) Mauna Loa, Hawaii, USA, (c) Tutuila, American Samoa, and (d) South Pole, Antarctica.

Figure
Figure3bshows the mean age of air of SF 6 between the northern and southern hemispheres and τ ex .The mean age of air is younger in MALTA (1.2 years) than in GEOS-Chem (1.7 years) (we use values from 1998 onward to allow a sufficient time lag).Also in Figure3bis the mean τ ex .MALTA has a faster exchange time (1.4 years) than GEOS-Chem (1.5 years), as was the case for the age of air of SF 6 metric.Yang et al. (2019) derive τ ex of around 1.4 years and a SF 6 of around 1.5 years from 21 3D transport models and τ ex and a SF 6 of around 1.4 years from observations.The age of air of SF 6 in MALTA is shorter than all models presented in Yang et al. (2019) but τ ex falls within the presented range of models.Yang et al. (2019) also derive these metrics for GEOS-Chem in their work using a finer resolution grid of 2.5°× 2°.The age of air of SF 6 for GEOS-Chem derived here falls broadly in line with that presented in their work, whereas τ ex derived inYang et al. (2019) is slightly shorter (around 1.4 years).The shorter timescale for a SF 6 in MALTA indicates that MATLA may have some shortcomings in simulating transport from northern hemisphere midlatitudes to the northern hemisphere tropics, that is, transport is too rapid(Yang et al., 2019).

Figure 2 .
Figure 2. The monthly mean mole fraction from (a) the zonal average of a 3D model (GEOS-Chem), (b) MALTA and (c) the 12-box model for January 2005 for the emissions scenario described in Section 4.1.

Figure 3 .
Figure 3.A comparison of the north to south transport timescale of SF 6 using GEOS-Chem and MALTA from 1995 to 2009.(a) The mean seasonal pattern of the interhemispheric exchange time, τ ex , of SF 6 .Shading shows the 10%-90% percentile of the variability for each month.(b) A comparison of τ ex and the age of SF 6 , where the squares show the mean monthly derived values over 1995-2009 and lines extend the 10%-90% percentile of the variability for all months.The age of air of SF 6 is the time lag for the mean mole fraction in the southern hemisphere to equal that in the northern hemisphere.The interhemispheric exchange, τ ex , is a similar metric considering the imbalance of emissions between hemispheres.

Figure 4 .
Figure 4.A comparison between the stratosphere-troposphere exchange of SF 6 from 1980 to 2008 in GEOS-Chem and MALTA using the same emissions.(a) The monthly mean stratosphere-troposphere exchange between 1980 and 2008, where the shading shows the 95% variability.(b) Box plots show the median, interquartile range and 1.5× the interquartile range of the stratosphere-troposphere exchange.

Figure 5 .
Figure 5.The global average partial pressure anomaly of CFC-11 in MALTA due to the influence of the QBO is shown by the shading.The magenta line shows the east-west equatorial wind anomaly at 20 hPa from observations above Singapore.

Figure 6 .
Figure 6.Emissions of CFC-11 from 2000 to 2019 estimated using MALTA and mole fractions simulated using TOMCAT and an inverse framework.The blue line shows the mean global emissions estimated using the MALTA with the one standard deviation uncertainty shown by the error bars.The dashed orange line show the prior mean emissions used in the emissions estimate with the one standard deviation uncertainty shown by the shading.The black line shows the annual total emissions used in TOMCAT to simulate the mole fractions used as input to the inverse framework.