Limits on the ability of global Eulerian models to resolve intercontinental transport of chemical plumes

Quasi-horizontal chemical plumes in the free troposphere can preserve their concentrated structure for over a week, enabling transport on intercontinental scales with important environmental impacts. Global Eulerian chemical transport models (CTMs) fail to preserve these plumes due to fast numerical dissipation. We examine the causes of this dissipation and how it can be cured. Goddard Earth Observing System (GEOS-5) meteorological data at 0.25× 0.3125 horizontal resolution and ∼ 0.5 km vertical resolution in the free troposphere are used to drive a worldwide ensemble of GEOS-Chem CTM plumes at resolutions from 0.25× 0.3125 to 4× 5, in both 2-D (horizontal) and 3-D. Two-dimensional simulations enable examination of the sensitivity of numerical dissipation to grid resolution. We show that plume decay is driven by flow divergence and shear, filamenting the plumes until GEOS-Chem’s highorder advection scheme cannot resolve gradients and fast numerical diffusion ensues. This divergence can be measured by the Lyapunov exponent (λ) of the flow. Dissipation of plumes is much faster at extratropical latitudes than in the tropics and this can be explained by stronger divergence. The plume decay constant (α) is linearly related to λ, and increasing grid resolution provides only modest benefits toward plume preservation. Three-dimensional simulations show near-complete dissipation of plumes within a few days, independent of horizontal grid resolution and even in the tropics. This is because vertical grid resolution is inadequate in all cases to properly resolve plume gradients. We suggest that finer vertical grid resolution in the free troposphere is essential for models to resolve intercontinental plumes, while current horizontal resolution in these models (∼ 1) is sufficient.

. Free tropospheric plumes resulting from stratospheric intrusions can retain 150 ppb of ozone over a period of weeks (Trickl et al. 2011). Such global-scale transport with little dilution has important implications for environmental impacts, interactions with weather, and chemical aging. 5 Eulerian models used for simulating global atmospheric transport fail to reproduce this persistent layered structure. The modeled plumes dissipate within days by mixing with the background (Heald et al. 2003;Vuolo et al. 2009). Eulerian models simulate transport as a flux divergence for fixed grid cells, and the rapid dissipation implies a large transportation error from numerical diffusion even when a highly accurate advection algorithm is used (Rastigejev et al. 2010). A Lagrangian approach, where transport is calculated for individual air parcels carried by with the flow with no interaction between neighbors, would 10 avoid this problem (Khosrawi et al. 2005). Global Lagrangian models have been used with success in the stratosphere to describe the sharp gradients at the edge of the polar vortex (Fairlie et al. 1999;Hoppe et al. 2014;Konopka et al. 2003).
However, they are not used in global tropospheric applications because of difficulties in dealing with convective motions, homogeneity of coverage, and non-linear chemistry (Brasseur and Jacob 2016). Adaptive mesh refinement techniques have shown promise in addressing this issue in Eulerian models (Semakin and Rastigejev, 2016) but are computationally complex . 15 There is a need to understand why persistent free tropospheric plumes are so rapidly dissipated in Eulerian models, and how this behavior can be corrected.
A theoretical study by Rastigejev et al. (2010) examined the causes of the fast numerical dissipation of intercontinental plumes in Eulerian models. Numerical diffusion is due to finite differencing of the advection equation on the model grid such that the 20 gradients between grid cells are imperfectly described. The order of a numerical advection scheme is defined by the number of adjacent grid cells used to resolve a local gradient. Rastigejev et al. (2010) showed that a highly accurate, third-order, finitevolume advection scheme such as is used in GEOS-Chem (Lin and Rood 1996) successfully preserves plume structures for over 10 days in a uniform flow, but fails rapidly when real-world divergence is applied. Flow divergence acts to filament, stretch, and thin the plume until it is resolved by only a few grid cells. At that point the gradient can no longer be represented 25 with a high-order scheme; the scheme collapses to first order, resulting in very fast numerical dissipation. Increasing grid resolution only delays the onset of this effect, and may amplify it through the additional flow divergence. Rastigejev et al. (2010) proposed that the plume dissipation rate in a stretched flow is eventually set by the Lyapunov exponent λ = ∂u/∂x of the flow, defined as the exponential rate at which adjacent trajectories (aligned with the wind speed vector u) diverge from each other. Increasing the grid resolution Δx of the model would then slow down the rate of decay only as Δx 0.5 , rather than 30 Δx 3 as might be expected from a third-order advection scheme. Given that computational costs in Eulerian models typically rise at a rate of Δx 2 to Δx 3 , this implies that a computationally expensive increase in resolution will yield only marginal reductions in the errors driving numerical plume dissipation.
In this paper, we examine whether the theory of Rastigejev et al. (2010) can explain the fast numerical decay of free tropospheric plumes in Eulerian models, and what the implications are for curing this problem through increasing grid resolution. We use for this purpose global 2-D (horizontal) and 3-D versions of GEOS-Chem to simulate atmospheric flow at horizontal grid resolutions ranging from 0.25°×0.3125° (~25×30 km 2 ) to 4°×5° (~400×500 km 2 ) and with a native vertical resolution of ~0.6 km in the free troposphere. We quantify the decay rate for plumes originating in different locations around 5 the world, relate this decay rate to flow stretching, and conclude as to the potential to preserve the plumes through the use of improved grids.

Theory
The theory presented by Rastigejev et al. (2010) for numerical diffusion of stretched plumes begins with the advection equation (1) 10 where n is the number density of an inert chemical ("tracer") and u is the wind vector. Rastigejev et al. (2010) expressed this equation in its advective form (2) and included a diffusivity term D to describe numerical diffusion: 15 Here C = n/na is the volume mixing ratio (VMR) and na is the number density of air. This form explicitly accounts for the effect that numerical diffusion has on the modeled flow. Without this numerical diffusion, the advection equation would strictly conserve the VMR. When numerical diffusion is included, the VMR decays over time as the plume dissipates.

20
Let us consider now the conceptual picture of a model plume with uniform VMR diluting by numerical diffusion into a background atmosphere with a VMR of zero. The plume has surface area S and volume V. Mass balance for the plume is given by where k is a unit vector normal to the plume surface. Rastigejev et al. (2010) defined a characteristic length rb for decay across the edge of the plume so that ∇C ~ C/rb. They further defined a characteristic width of the plume as W = V/S. Thus equation This implies an exponential decay in C, such that where the decay constant ɑ is given by 10 = and Δt is some time interval.
The decay rate of the plume is proportional to the numerical diffusivity D, which is dictated by the order f of accuracy of the numerical advection scheme: D ~ Δx f . However, the decay rate also depends on rb. In a divergent flow, stretching of an initially 15 broad plume (W >> rb) causes rb to decrease, and hence the decay rate to increase. This stretching can be represented by the Lyapunov exponent, defined as where (u,v) are the (x,y)-components of the velocity and x is taken as the direction of stretching. The Lyapunov exponent 20 defines the exponential rate constant at which initially adjacent trajectories diverge.
Stretching in a divergent flow thins the plume, while numerical diffusion thickens it. Under constant divergence (constant λ), an equilibrium size for rb is reached when these two processes proceed at the same rate. Since the rate constant for diffusion is ~D/rb 2 , and the rate constant for stretching is λ, equilibrium is reached when 25 Atmos. Chem. Phys. Discuss., doi:10.5194/acp-2016-943, 2016 Manuscript under review for journal Atmos. Chem. Phys. Published: 24 November 2016 c Author(s) 2016. CC BY 3.0 License.

=
Replacing into (6), we find Thus the rate of decay in a stretched flow is less sensitive to D than expected. Furthermore, if the plume has stretched to be only a few gridboxes thick, then the gradient across the plume boundary cannot be resolved by a high-order scheme anymore and any numerical advection scheme collapses to first order (Godunov 1959). Under these conditions D ~ Δx and thus ɑ ~ Δx 0.5 ; the decay rate improves only as the square root of the grid resolution.

10
We also see from equation (9) that the decay rate increases with the rate of stretching as measured by λ. Regions of divergent flow are expected to experience faster plume decay. Eventually, for a fully stretched plume we have W = rb. Under these conditions, equations (8) and (9) yield ɑ = λ and the decay rate is independent of the grid resolution -a remarkable result.
The Rastigejev et al. (2010) theory thus paints the following picture for the model decay of a free tropospheric plume in a 15 divergent flow (as is realistically found in the atmosphere) and its dependence on grid resolution Δx. A plume that is initially well-resolved on the model grid will decay with a rate constant propotional to Δx f , where f is the order of accuracy of the numerical advection scheme. As the plume stretches, shears and filaments, it becomes poorly resolved on the model grid, and at that point the decay rate becomes proportional to (λΔx) 0.5 , i.e. only weakly responsive to increasing grid resolution and dependent on the stretching rate. In fact, increasing grid resolution may increase the stretching of the flow by introducing 20 additional convergence-divergence zones that would be averaged out at coarser resolution. Under those conditions Rastigejev et al. (2010) find that the decay rate may be proportional to Δx 0.25 , an even weaker grid resolution dependence. Eventually, the filamented plume decays with a rate constant defined by λ and at that point very fast dissipation takes place that is resolutionindependent. This theory, if correct, has major implications for understanding the decay of free tropospheric plumes in Eulerian models, and the value of increasing grid resolution. In what follows we test the theory using global simulations in actual 25 atmospheric flow with the GEOS-Chem CTM.

Testing theory with the GEOS-Chem CTM
We simulate transport of free tropospheric plumes in v11-01e of the global Eulerian GEOS-Chem CTM originally described by Bey et al. (2001). The model is driven by winds and other meteorological data archived every 3 hours from the Data Assimilation System of the NASA Goddard Earth Observing System (GEOS-5) with 0.25°×0.3125° horizontal resolution on 72 vertical levels. The vertical grid resolution in the free troposphere between 4 and 8 km altitude is about 0.6 km. We apply 5 the model to an inert chemical tracer with only advection enabled. Subgrid transport processes, including convection and boundary-layer mixing, are disabled. Thus the model only solves for advection (equation (1)), using the 3-hour GEOS-5 FP archive of mean horizontal winds and instantaneous surface pressure. Horizontal winds are adjusted with a "pressure fixer" (Horowitz 2003) to ensure consistency with the 3-hour pressure change. Vertical winds are derived from divergence of the horizontal winds and the change in surface pressure. 10 Horizontal advection is calculated using the FFSL-3 finite volume scheme developed by Lin and Rood (Lin and Rood 1996) and commonly called "tpcore". This scheme uses the monotonic piecewiseparabolic method (PPM) when the Courant-Friedrichs-Lewy number (CFL) is less than or equal to one, and a semi-Lagrangian method for CFL > 1. A semimonotonic PPM is used in the vertical direction with the enforcement of Hyunh's second monotonicity constraint. The 15 FFSL-3 scheme is formally third-order accurate in space, such that increasing the grid resolution Δx by a factor of 2 should reduce numerical errors by a factor of 8.
We conduct 2-D (horizontal) and 3-D simulations at 5 different horizontal resolutions: 0.25°×0.3125° (native), 0.5°×0.625°, 1°×1.25°, 2°×2.5°, and 4°×5°. 2-D simulations allow an analysis of the effect of horizontal resolution over a factor of 16 (from 20 0.25°×0.3125° to 4°×5°) to test the theoretical dependences of plume dissipation on grid resolution derived by Rastigejev et al. (2010). We cannot carry out a similar analysis in the vertical because the ~0.6 km native vertical resolution of the GEOS-5 data in the free troposphere is too coarse. In all cases, winds and pressures are retrieved from the GEOS-5 FP archive at the 0.25°×0.3125° native resolution and subsequently averaged spatially to drive coarser-resolution simulations as is routinely done in GEOS-Chem (Bey et al. 2001;Philip et al. 2015). A dynamical timestep of 5, 5, 10, 15, and 30 minutes respectively 25 is used for each resolution. Decreasing the timestep to 5 minutes for all simulations had no significant effect. All concentrations shown and discussed are based on instantaneous VMRs stored at one-hour intervals.
2-D simulations are performed by taking the pressure-weighted average of the wind velocity in each atmospheric column and setting the surface pressure tendency to zero. Although clearly idealized, there is some realism to the 2-D simulations in that 30 free tropospheric layers are vertically stratified and most of the shearing and dissipation can be expected to take place in the horizontal. Most relevantly, the 2-D simulations allow us to test the theory of Rastigejev et al. (2010) for the sensitivity of plume dissipation to grid resolution.
We conduct simulations of the first 9 days of July 2015 for plumes initialized in different locations around the world with a homogeneous unit mixing ratio over a cuboid 12° in latitude by 15° in longitude, and zero outside. This size is chosen so that the initial plume is coarsely resolved at 4°×5° but finely resolved at 0.25°×0.3125°. 90 non-interacting plumes are initialized over the global domain ( Fig. 1) to examine how location affects numerical diffusion. In our 2-D simulations, only horizontal 5 advection is enabled, and the model domain is restricted to a single vertical layer. In our 3-D simulations, vertical advection is enabled, and each plume is initialized at 3.9 km pressure altitude over a single GEOS-5 model level that is 470 m thick.

Numerical diffusion
Exact solution of the advection equation translates mixing ratios downwind without altering them. In other words, initial mixing ratios in a plume remain constant as the plume is advected downwind. Any plume decay in our advection-only simulation must be the result of numerical diffusion. In the real atmosphere, plumes decay by molecular diffusion that operates 15 on millimeter scales and is the end result of the turbulent eddy cascade that filaments the plume into finer and finer strands (Brasseur and Jacob 2016). This subgrid turbulence is particularly fast for vertical mixing in the boundary layer, and is typically parameterized in models with a turbulent diffusion scheme. It is usually ignored in the horizontal direction or in the free troposphere, under the assumption that spurious numerical diffusion effectively carries out the mixing. Numerical diffusion arises from finite differencing over grid cells when solving the advection equation. Odd-order schemes such as the PPM tend to introduce diffusion, artificially smoothing the solution in areas with sharp concentration gradients.
Even-order schemes tend instead to be dispersive, producing spurious oscillations, and this artifact information is even less desirable than numerical diffusion. Higher-order schemes also tend to produce spurious oscillations when there are discontinuities in the concentration field (Godunov 1959;Brasseur and Jacob 2016). To get around this, modern advection 5 schemes such as FFSL-3 employ flux limiters that locally reduce the scheme to first order in the vicinity of discontinuities.
This prevents spurious oscillations at the cost of increasing numerical diffusion.
Numerical diffusion in GEOS-Chem is illustrated in Fig. 2 with an example of a 2-D plume at 1°×1.25° grid resolution. The plume decays with time due to numerical diffusion. This decay is reflected by an increase in the plume area and a decrease in 10 the plume VMR. The rate of decay can be measured by the reduction in the plume maximum VMR, as done by Rastigejev et al. (2010) and expressed in equation (5) in terms of an exponential plume decay constant ɑ. In situations where the plume is resolved by a large number of grid cells, the maximum VMR can be buffered from the effects of numerical diffusion by surrounding grid cells, even as the plume frays at its edges. In the example shown here, the maximum value remains nearly unchanged for 4 days as a result of this buffering and hence ɑ is near zero. Even beyond 4 days, ɑ can be highly variable 15 depending on the local flow divergence (Rastigejev et al., 2010) and tends to decrease as the plume dissipates because dissipation smooths the VMR gradient.
An alternate metric of numerical diffusion is the size of the plume. The thick red contour in Fig. 2 shows the minimum area containing 90% of the total mass of tracer in the plume. As the simulation progresses, diffusion of the plume increases this 20 area. We define the square root of this area as the characteristic size of the plume, normalized by the value at plume initialization. In 3-D, the plume size is taken as the cube root of the total volume occupied by 90% of the tracer mass, after accounting for differences in air density. Plume size is a more sensitive indicator of plume diffusion, as shown for our example in the bottom panel of Fig. 2, because it accounts for the fraying at the plume edges and is a smoother function than ɑ. Here we will mostly use the decay constant ɑ of maximum VMR as metric of plume decay, for consistency with Rastigejev and to 25 compare to theory (Section 2), but we will also show some results for plume size.
Atmos. Chem. Phys. Discuss., doi:10.5194/acp-2016-943, 2016 Manuscript under review for journal Atmos. Chem. Phys. Published: 24 November 2016 c Author(s) 2016. CC BY 3.0 License. Figure 2. Numerical diffusion of a 2-D (horizontal), inert plume in GEOS-Chem for a 9-day simulation at 1°×1.25° horizontal resolution. This particular plume was initialized between Australia and New Zealand (Fig. 1). The red contour shows the minimum area containing 90% of the tracer mass. From top to bottom, the lower three panels show: the normalized maximum VMR in the plume; the 6-hour moving average decay constant ɑ calculated from equation (5); and the plume size, defined as the square root of 5 the 90% contour area and normalized by the initial value.

Stretching
Plume stretching can be quantified by the local Lyapunov exponent of the flow, as defined in equation (6) Rearranging equation 11 gives This can be directly applied in an Eulerian model framework, acknowledging the separate treatment of winds in each 5 dimension. For taxicab geometry, as in the orthogonal latitude-longitude discretization, we measure the separation δ(t) as (Δx + Δy) at time t = 0, where Δ denotes the grid cell spacing. The separation at time t + Δt is then given by

| ( + Δ )| = (Δ + |Δ Δ |) + (Δ + |Δ Δ |)
Here, Δu and Δv refer to the change in wind speed between the current grid cell and the downstream cell in each direction. 10 Using the absolute value prevents flow stretching in one direction from being offset by compression in another. Instead, this metric responds to flow stretching in either the x-or y-direction. We also assume that the change in separation will be small relative to the initial grid spacing, allowing us to approximate ln(x) ≈ (1-x). Replacing into equation (12) At each grid cell, the above equation is applied to yield an estimator for the rate of horizontal flow stretching. Figure 1 displays the mean Lyapunov exponents at 0.25°×0.3125° for the full 9-day simulation period in the 2-D flow. As discussed by Stohl (2001), the weaker synoptic-scale eddies at low latitudes result in less flow stretching relative to the higher latitudes. The rate of plume decay increases with latitude. A tropical plume with initial area of 12°×15° retains over 99% of its original maximum VMR after 9 days at a resolution of 0.25°×0.3125°, 98% at 0.5°×0.625°, and 89% at 1°×1.25°. Outside of the tropics, these values fall to 82% at 0.25°×0.3125°, 59% at 0.5°×0.625° and 38% at 1°×1.25°. A plume in the tropics is better preserved at 1°×1.25° than a plume outside the tropics at 0.25°×0.3125°. 5 Rastigejev et al. (2010) presented a single example of a Chinese plume transported over the Pacific in 2-D flow as illustration of their theory. Starting from the same initial 12°×15° plume dimension, they found the maximum VMR to drop to 10% of its original value after 9 days 1°×1.25°. Our results considering a large ensemble of plumes do not show such a drastic decay. In fact, it would seem that the 0.25°×0.3125° resolution is largely successful at preserving plumes over the 9-day period. As we 10 will see in section 6, this success does not hold for 3-D plumes; but we focus on 2-D plumes in this section to better understand dependences on flow stretching and grid resolution.

Relationship to flow divergence
Following the theory of Rastigejev et al. (2010) as summarized in Section 2, we examined the relationship between plume decay (as represented by the plume decay constant ɑ) and flow stretching (as represented by the Lyapunov exponent λ).  Fig. 4 show in general a strong correlation between ɑ and λ, demonstrating that flow stretching plays a major role in driving plume decay. The relationship is linear at 4°×5° grid resolution, as would be expected for a fully stretched plume (section 2). At the higher resolutions, ɑ remains near zero when λ is low (tropical plumes) implying that the plume remains well-defined when stretching is weak. However, when stretching is strong (λ > 10 -5 s -1 ), we still find a linear relationship between ɑ and λ, indicating fast plume decay from plume filamentation even at the higher resolutions. Rastigejev et al. (2010) 10 gave ɑ = λ for the fully stretched plume on the model grid, but we find weaker slopes that decrease with increasing resolution (Fig. 4). This suggests that the model actually retains some capability to describe cross-plume gradients, and the Rastigejev et al. (2010) assumption of a sharp discontinuity on the model grid must be viewed as a limiting case.

Results in
The difference in the slope of the regression lines gives an approximate measure of the improvement gained by increasing 15 resolution by a factor of 4. Increasing resolution provides a "delay" in terms of the minimum λ at which the plumes being to decay rapidly due to stretched-flow diffusion. At high resolution and low λ, the maximum VMR is preserved for the full 200 hours, and minimal numerical diffusion occurs. At high λ or low resolution, the buffering is insufficient, and numerical diffusion proceeds at a rate that decreases by about a factor of 2 for every 4-fold increase in resolution This supports a central result of Rastigejev's theory that the plume decay rate decreases as the square root of the grid resolution. 20 The right-hand panel of Fig. 4 shows the response of the plume size to diffusion after 200 hours. At 4°×5°, the average size increase after 200 hours is a factor of 4.5, compared to 2.0 at 0.25°×0.3125°. Unlike the decay constant ɑ which relates to the maximum VMR, the plume size consistently increases as λ increases, reflecting the effect of stretching in fraying the plume edges even if the plume maximum VMR is preserved. The rate of improvement in plume size with resolution is therefore 25 slower and there is no minimum value of λ, as numerical diffusion affects the plumes immediately at all resolutions. Overall, a factor of 16 increase in grid resolution yields a reduction in plume size by a factor of about 3, compared to a factor of 4 for the decay constant.
Atmos. Chem. Phys. Discuss., doi:10.5194/acp-2016-943, 2016 Manuscript under review for journal Atmos. Chem. Phys.   (Fig. 1). Each datapoint corresponds to a single plume and the grid resolution is identified by color. Linear regression lines are calculated using reduced major axis regression, discounting points below the "regression cutoff"

Effective order of accuracy
The PPM advection scheme used in GEOS-Chem is 3 rd -order accurate, meaning that the accuracy should improve as Δx 3 for well-resolved plumes. As we have seen above, this is far from the case for a long-lived plume in a divergent flow. The theoretical analysis of Rastigejev et al. (2010) indicates that the decay rate of a well-resolved plume should initially improve 10 with the order of accuracy of the advection scheme, but that the rate of improvement should fall off to Δx 0.25-0.5 as plume stretching limits the ability to resolve gradients. Furthermore the decay rate is expected to eventually become grid-independent as the dimension of the filamented plume becomes comparable to the model grid (section 2). Here we use our ensemble of 2-D simulations ranging over a factor of 16 in grid resolution to evaluate that result, taking the plume decay rate constant ɑ as a measure of accuracy. 15 Figure 5 shows the average improvement (reduction) in ɑ with increasing grid resolution, for the plumes in different latitude bands and as a function of plume age. In the tropics, there is in general little flow divergence and numerical diffusion is weak as a result. Figure 5 shows that ɑ in the fresh plume improves as Δx 3 , consistent with the third-order accuracy of the PPM Atmos. Chem. Phys. Discuss., doi:10.5194/acp-2016-943, 2016 Manuscript under review for journal Atmos. Chem. Phys. Published: 24 November 2016 c Author(s) 2016. CC BY 3.0 License. advection scheme, and even in the aged plume the improvement scales as Δx 2 . Thus we find that flow divergence does not limit the gains from increasing grid resolution, at least for these 2-D plumes. As shown in Fig. 3, a 1°×1.25° grid resolution seems sufficient to simulate long-range transport of tropical plumes with little numerical diffusion.
Outside the tropics, where flow divergence is greater, the effective order of accuracy is smaller and shows greater variability 5 between grid resolutions with plume age. It starts second-order (Δx 2 ) but decreases as the plume ages and filaments. By day 5-6, the average order of convergence has decreased to Δx 0.5 for grid resolutions coarser than 1°×1.25°. At higher resolutions the rate of improvement with resolution is greater, albeit still of the order of Δx. There is curvature in the response of the plume decay constant to grid resolution, such that the benefit of increasing grid resolution increases as the resolution gets finer. This is again in agreement with the theory of Rastigejev et al. (2010), where the purpose of increasing grid resolution is to drop 10 below the scale where plume gradients are well resolved.  Atmos. Chem. Phys. Discuss., doi:10.5194/acp-2016-943, 2016 Manuscript under review for journal Atmos. Chem. Phys. Published: 24 November 2016 c Author(s) 2016. CC BY 3.0 License.

3-D plume decay
We now turn to 3-D simulations for a more practical evaluation of the gains that could be made from increasing model resolution. An important distinction here is that we cannot explore a wide range of grid resolutions in the vertical; the ~0.6 km native vertical resolution of the GEOS-5 data in the free troposphere is comparable to the observed thickness of free tropospheric plumes (Thouret et al. 2000). Such coarse resolution in the free troposphere seems typical of the current generation 5 of models, which have emphasized improving horizontal resolution more than vertical resolution. For example, the ERA-Interim re-analysis, produced by the European Center for Medium-range Weather Forecasts (ECMWF), has a similar mean vertical resolution of 570 m in the free troposphere (Dee et al. 2011).
Global mean plume decay rates and plume sizes in 3-D are shown in Fig. 6 as a function of plume aging times, and compared 10 to the 2-D cases discussed previously. In the 2-D simulations, each doubling of resolution yields a 10-20% improvement in the final maximum VMR after 9 days, up to a value of 89% at a resolution of 0.25°×0.3125°. The plume size improvement with increasing resolution is smaller but equally consistent. After 9 days, the plume size is double its initial size at 0.25°×0.3125°. In 3-D the numerical diffusion is considerably larger. At the 0.25°×0.3125° grid resolution, the maximum VMR after 9 days drops to 13% of its original value, and the plume size increases by a factor of 6 from its original value. The 15 reason is that the ability to preserve the plume is limited by numerical diffusion in the vertical, where the plume is initially poorly resolved in all cases because of the low native vertical resolution in the GEOS-5 fields.
There is also a counterproductive aspect to increasing horizontal resolution in a 3-D simulation. As the horizontal resolution is increased, fine-scale vertical eddies are resolved that increase the vertical stretching of the plume, compromising the 20 advantages gained from the slower horizontal diffusion. This is highlighted by the negligible improvement in plume size after 9 days between 3-D simulations at grid resolutions of 0.5°×0.625° and 0.25°×0.3125°. The increased vertical diffusion almost completely offsets the improvement provided by reducing spurious horizontal diffusion.
Atmos. Chem. Phys. Discuss., doi:10.5194/acp-2016-943, 2016 Manuscript under review for journal Atmos. Chem. Phys. Published: 24 November 2016 c Author(s) 2016. CC BY 3.0 License. In 2-D, increasing resolution yields consistent benefits, such that simulations at or finer than 1°×1.25° achieve near-zero diffusion in the tropics. In 3-D, high rates of flow stretching remain correlated with higher rates of diffusion, but increasing 10 the horizontal resolution yields a smaller relative decrease in the rate of diffusion, and the overall rate of convergence is of the order of Δx 0.25 or worse. We find that extratropical regions still consistently experience greater rates of numerical diffusion.
However, even in the tropics where flow stretching is slow and 2-D simulations performed well, vertical diffusion due to poor vertical resolution provides a lower bound of 3×10 -6 s -1 for ɑ, corresponding to a decay time scale of 3 days. This shows that the ability of current global Eulerian models to resolve free tropospheric plumes is limited by the vertical grid resolution. 15 Assuming isotropy in requirements between the horizontal resolution (which we were able to investigate in detail) and the vertical resolution, we conclude that a 100-m vertical resolution would be required to preserve the transport of free tropospheric plumes.

7 Conclusions
We examined why global models are unable to simulate the intercontinental-scale transport of quasi-horizontal plumes in the free troposphere, dissipating them in a few days instead of preserving their coherence. Our focus was to test theoretical results by Rastigejev et al. (2010) that this dissipation is due to fast numerical diffusion in the divergent flow typical of the free troposphere. Divergence (shear, stretching) causes the plume to filament rapidly to the point when it is not properly resolved 10 on the model grid. At that point fast dissipation takes place in the model regardless of the accuracy of the advection scheme and only weakly dependent on grid resolution.
We conducted a large worldwide ensemble of simulations of free tropospheric plumes with the GEOS-Chem chemical transport model driven by NASA GEOS-5 assimilated meteorological data. The simulations used horizontal resolutions 15 ranging from 0.25°×0.3125° (native GEOS-5) to 4°×5°, including only advection (no subgrid turbulence or chemistry).
Restriction to advection allowed us to focus on numerical diffusionin a purely advective problem the plume should not dilute, even in a divergent flow. We diagnosed plume decay caused by numerical diffusion in the model by the decrease in maximum volume mixing ratio (VMR) in the plume (decay rate constant α) and by the increase in plume size. Native vertical resolution in GEOS-5 (and other current meteorological models) is ~0.5 km in the free troposphere, too coarse to adequately 20 resolve vertical gradients in plumes (typically ~1 km thick). Thus we conducted both 2-D (horizontal) and 3-D simulations.
Restriction to 2-D allowed us to investigate in detail the sensitivity to grid resolution for initially well-resolved plumes (12°×15°). Extension to 3-D allowed us to examine numerical diffusion in realistic model situations.
We find that extratropical plumes decay much faster than tropical plumes, and that this can be explained by stronger flow divergence measured by the Lyapunov exponent of the flow (λ). Under strongly divergent flow typical of the extratropics, the 5 rate of plume decay varies linearly with λ and improves only as the square root of the grid resolution Δx. We find that the sensitivity of the plume decay to grid resolution decreases as the plume ages, initially improving as Δx 3 (the order of accuracy of the GEOS-Chem advection scheme), and eventually decaying after a few days to Δx 2 in the tropics and Δx 0-1 in the extratropics.
10 3-D plume decay in our simulations is much faster than in 2-D, and consistent with the general inability of models to preserve the coherence of free tropospheric plumes. The plume decay rate in 3-D still depends on horizontal flow divergence, but the sensitivity to horizontal grid resolution is weaker and the decay is instead limited by the coarse vertical resolution. Vertical numerical diffusion is very fast, and is amplified at finer horizontal resolution by vertical eddies that would be smoothed out at coarser horizontal resolution. Even tropical plumes decay with a time constant of about 3 days. 15 Thus we find that increasing vertical grid resolution in the free troposphere to ~100 m is an essential first step for models to resolve the intercontinental-scale transport of free tropospheric plumes. Increasing horizontal resolution beyond 1° is futile.
Even then, the modeling problem will remain challenging in strongly divergent flow typical of high latitudes. More advanced solutions might involve adaptive grids designed to resolve local conditions of large chemical gradients and large divergence, 20 or embedding Lagrangian plumes in the global Eulerian modeling framework.

Author contributions
SDE and DJJ designed the experiments. SDE developed model code and performed all experiments. SDE prepared the manuscript in collaboration with DJJ.

Competing interests 25
The authors declare that they have no conflict of interest.