Development and technical paper 30 Apr 2021
Development and technical paper  30 Apr 2021
Using the Després and Lagoutière (1999) antidiffusive transport scheme: a promising and novel method against excessive vertical diffusion in chemistrytransport models
 ^{1}LMD/IPSL, École Polytechnique, Institut Polytechnique de Paris, ENS, PSL Research University, Sorbonne Université, CNRS, Palaiseau France
 ^{2}École des PontsParisTech, MarnelaVallée, France
 ^{a}currently at: Aria Technologies, BoulogneBillancourt, France
 ^{1}LMD/IPSL, École Polytechnique, Institut Polytechnique de Paris, ENS, PSL Research University, Sorbonne Université, CNRS, Palaiseau France
 ^{2}École des PontsParisTech, MarnelaVallée, France
 ^{a}currently at: Aria Technologies, BoulogneBillancourt, France
Correspondence: Sylvain Mailler (sylvain.mailler@lmd.polytechnique.fr)
Hide author detailsCorrespondence: Sylvain Mailler (sylvain.mailler@lmd.polytechnique.fr)
The potential of the antidiffusive transport scheme proposed by Després and Lagoutière (1999) for resolving vertical transport in chemistrytransport models is investigated in an idealized framework with very encouraging results. We show that, compared to classical higherorder schemes, the Després and Lagoutière (1999) scheme reduces numerical diffusion and improves accuracy in idealized cases that are typical of atmospheric transport of tracers in chemistrytransport models. The increase in accuracy and the reduction in diffusion are substantial when, and only when, vertical resolution is insufficient to properly resolve vertical gradients, which is very frequent in chemistrytransport models. Therefore, we think that this scheme is an extremely promising solution for reducing numerical diffusion in chemistrytransport models.
Reducing numerical errors in chemistrytransport models (CTMs) is a necessary task for those wishing to improve these models' performance. Among the wellknown errors in Eulerian CTMs, excessive numerical diffusion in all directions is a wellknown drawback of many such models, and in the past decade, excessive or poorly represented vertical transport and diffusion has been identified as a major cause for numerical dispersion in these models (Vuolo et al., 2009; Emery et al., 2011; Mailler et al., 2017). This excessive vertical diffusion may have a strong impact on the representation of groundlevel ozone concentrations due to spurious transport of stratospheric ozone into the troposphere (Emery et al., 2011) and also hinders the ability of Eulerian CTMs to represent accurately intercontinental transport of densely polluted plumes such as volcanic plumes (Colette et al., 2010; Mailler et al., 2017; Lachatre et al., 2020). While CTMs manage to represent such plumes in terms of general location, they typically fail to maintain the finescale structure of the plumes and tend to dilute them too much compared to observations (Colette et al., 2010; Mailler et al., 2017).
Among other more modelspecific solutions, Emery et al. (2011) suggest as the main paths to solve this problem increasing vertical resolution and improving the vertical advection scheme, in their case switching from a firstorder to a secondorder advection scheme. In the same line, Eastham and Jacob (2017) and Zhuang et al. (2018) have discussed the need for increased vertical resolution in order to adequately represent longrange advection of chemical plumes.
In the line of improving the vertical advection scheme, Lachatre et al. (2020) describe the implementation of the Després and Lagoutière (1999) antidiffusive transport scheme for vertical transport in the CHIMERE chemistrytransport model (Mailler et al., 2017) and its application to modeling the 18 March 2012 eruption of Mount Etna (Italy). They show that using this transport scheme reduces numerical diffusion and permits a better representation of the volcanic plume after longrange advection, thereby showing that it is possible to strongly reduce numerical diffusion in CTMs without increasing the number of vertical levels. However, that study, set in the framework of a fully fledged chemistrytransport model fed by reallife atmospheric fields, makes it difficult to fully disentangle the effect of the transport scheme itself from other effects such as uncertainties in emission fluxes and mass–wind inconsistencies in the forcing meteorological fields.
Therefore, the present study aims at answering the questions on the use of Després and Lagoutière (1999) that could not be addressed in the realistic framework of Lachatre et al. (2020). For that purpose, we have designed three idealized test cases that permit comparison of the performance of the Després and Lagoutière (1999) scheme with the classical schemes of Van Leer (1977) and Colella and Woodward (1984), not only in terms of diffusion but also in terms of accuracy compared to the exact solution, which could not be done in Lachatre et al. (2020). We also examine how the performance of Després and Lagoutière (1999) compares to the two abovecited schemes and to the order1 upwind Godunov scheme (Godunov and Bohachevsky, 1959) when resolution increases.
Section 2 describes the numerical methods that have been used, their implementation, the discretization strategies and the cases that have been designed for the study. Section 3 presents simulation outputs and the diagnostics that have been designed to compare the different transport strategies with each other. In Sect. 4 the results are discussed, and Sect. 5 shows our conclusions.
Continuity equation for the motion of air is as follows:
where u represents wind speed, C is air concentration (all species together) in molecules per unit volume and Φ=Cu represents the air flux vector.
The continuity equation for species s is as follows:
where C_{s} is concentration of species s in molecules per unit volume and Φ_{s}=C_{s}u is the flux vector for species s. Equivalently, Eq. (2) becomes
where α_{s} is the mixing ratio for species s:
Chemistrytransport models try to solve Eq. (3) as accurately as possible on their discretized grids, while keeping the cost of numerical resolution under control.
2.1 1D discretization of advection and advection schemes
In 1D, Eq. (2) becomes
Here we follow a semiLagrangian approach by identifying the air parcels that have entered cell i through its left boundary, and conversely the air parcels that have left cell i through its right boundary (we suppose that the wind is positive). Let Δ^{−}x be so that the Lagrangian trajectory starting at ${x}_{i\frac{\mathrm{1}}{\mathrm{2}}}{\mathrm{\Delta}}^{}x$ at time t passes through x_{i−½} at time t+Δt. Δ^{−}x is the distance traveled by the last air particle entering grid cell i at time t+Δt. Let us define the wind speed representative of facet $i\frac{\mathrm{1}}{\mathrm{2}}$ between t and t+Δt as ${\stackrel{\mathrm{\u203e}}{u}}_{i\frac{\mathrm{1}}{\mathrm{2}}}=\frac{{\mathrm{\Delta}}^{}x}{\mathrm{\Delta}t}$. If we define the wind speed representative of right facet ${\stackrel{\mathrm{\u203e}}{u}}_{i+\frac{\mathrm{1}}{\mathrm{2}}}$ in a similar way to $\frac{{\mathrm{\Delta}}^{+}x}{\mathrm{\Delta}t}$ where Δ^{+}x is so that the Lagrangian trajectory starting at ${x}_{i+\frac{\mathrm{1}}{\mathrm{2}}}{\mathrm{\Delta}}^{+}x$ at time t passes through x_{i+½} at time t+Δt, then Eq. (5) can be discretized over cell i as follows:
where
and
Equation (6) is verified exactly with no particular hypothesis on the wind speed u(x,t) nor the concentration field C_{s}(x,t). For air concentration, the continuity equation can be discretized in the same way:
In order to permit monotonicity of the advection scheme in terms of mixing ratios (i.e., the mixing ratio for species s stays within its initial range), we reformulate Eq. (6) by using fluxes and mixing ratios instead of winds and concentrations by introducing
and
With these notations, Eqs. (6) and (9) become
and
Equations (12) and (13) are a fluxform reformulation of semiLagrangian Eqs. (6)–(9). The form of Eqs. (12) and (13) makes it straightforward to verify that if Eq. (13) is verified, and if ${\stackrel{\mathrm{\u203e}}{\mathit{\alpha}}}_{s,i\frac{\mathrm{1}}{\mathrm{2}}}$ lies between ${\stackrel{\mathrm{\u203e}}{\mathit{\alpha}}}_{s,i\mathrm{1}}$ and ${\stackrel{\mathrm{\u203e}}{\mathit{\alpha}}}_{s,i}$ (and if ${\stackrel{\mathrm{\u203e}}{\mathit{\alpha}}}_{s,i+\frac{\mathrm{1}}{\mathrm{2}}}$ lies between ${\stackrel{\mathrm{\u203e}}{\mathit{\alpha}}}_{s,i}$ and ${\stackrel{\mathrm{\u203e}}{\mathit{\alpha}}}_{s,i+\mathrm{1}}$) then the resulting advection scheme guarantees monotonicity of mixing ratios, which is physically desirable and is the reason why chemistrytransport models usually resolve advection of trace species using an approach based on Eq. (12) rather than a straightforward resolution of Eq. (6).
Regarding chemistrytransport models, in practice, approximated values of ${\stackrel{\mathrm{\u203e}}{F}}_{i\frac{\mathrm{1}}{\mathrm{2}}}$ are inferred from the wind and density fields provided by the forcing meteorological dataset. Here we will avoid the typical problem of masswind inconsistencies discussed in, for example, Jöckel et al. (2001), Emery et al. (2011) and Lachatre et al. (2020), by working with analytically defined nondivergent mass fluxes, and constant and uniform air density, so that Eq. (12) is verified exactly by construction.
This simplified framework will permit us to focus on the transport scheme of the chemistrytransport model whose task, in fluxform, is to estimate the values of ${\stackrel{\mathrm{\u203e}}{\mathit{\alpha}}}_{s,i\pm \frac{\mathrm{1}}{\mathrm{2}}}$ that are needed for numerical resolution of Eq. (12).
2.1.1 Advection schemes and tracer flux calculation
The Godunov donorcell scheme
The most simple way of estimating ${\stackrel{\mathrm{\u203e}}{\mathit{\alpha}}}_{s,i\pm \frac{\mathrm{1}}{\mathrm{2}}}$ is the Godunov donorcell scheme (adapted from Godunov and Bohachevsky, 1959), simply evaluating ${\stackrel{\mathrm{\u203e}}{\mathit{\alpha}}}_{s,i,k+\frac{\mathrm{1}}{\mathrm{2}}}$ as follows:
This order1 scheme is cheap, robust, linear, monotonous and massconservative but extremely diffusive. It is therefore important to find more accurate ways to estimate ${\stackrel{\mathrm{\u203e}}{\mathit{\alpha}}}_{s,i+\frac{\mathrm{1}}{\mathrm{2}}}$.
The Van Leer (1977) scheme
The secondorder slopelimited scheme of Van Leer (1977) brought to our notations and assuming uniform air concentration yields the following expression of ${\stackrel{\mathrm{\u203e}}{\mathit{\alpha}}}_{s,i+\frac{\mathrm{1}}{\mathrm{2}}}$ (for ${\stackrel{\mathrm{\u203e}}{F}}_{i+\frac{\mathrm{1}}{\mathrm{2}}}>\mathrm{0}$):
where $\mathit{\nu}=\frac{{\mathrm{\Delta}}^{+}x}{{x}_{i+\frac{\mathrm{1}}{\mathrm{2}}}{x}_{i\frac{\mathrm{1}}{\mathrm{2}}}}$ is the Courant number for the donor cell. If ν>1, then more mass leaves the cell than the mass that was initially present and the Courant–Friedrichs–Lewy condition is violated, yielding numerical instability. If α_{s,i} is a local extremum of mixing ratio ($({\mathit{\alpha}}_{s,k}{\mathit{\alpha}}_{s,k\mathrm{1}})({\mathit{\alpha}}_{s,k+\mathrm{1}}{\mathit{\alpha}}_{s,k})\le \mathrm{0}$), no interpolation is performed and ${\stackrel{\mathrm{\u203e}}{\mathit{\alpha}}}_{s,i+\frac{\mathrm{1}}{\mathrm{2}}}={\mathit{\alpha}}_{s,i}$ is imposed: in this case, the scheme falls back to the simple Godunov donorcell formulation (Eq. 14). This order2 scheme has been used for decades in chemistrytransport modeling, being a good tradeoff between reasonably weak diffusion, at least compared to more simple schemes such as the Godunov donorcell scheme, and small computational burden compared to higherorder schemes such as the piecewise parabolic method (Colella and Woodward, 1984).
The Colella and Woodward (1984) piecewise parabolic method
The Colella and Woodward (1984) piecewise parabolic method (PPM) consists of performing a parabolic reconstruction of the concentration field inside each model cell using information from three upwind cells and two downwind cells, and applying limiters to preserve the scheme's monotonicity and stability. The detailed procedure is described in the seminal Colella and Woodward (1984) paper. Our implementation of this method has been adapted from the CASTRO Compressible Astrophysical Solver Almgren et al. (2010). While thirdorder by design, application of limiters in the vicinity of extrema introduces firstorder truncation errors in their vicinity so that thirdorder convergence is not expected with the PPM scheme as described in Colella and Woodward (1984) (Colella and Sekora, 2008). However, this limitation does not prevent the Colella and Woodward (1984) method giving much better results than simpler order2 schemes such as Van Leer (1977), so that the Colella and Woodward (1984) PPM scheme has been used successfully for a wide range of applications including meteorological modeling (Carpenter et al., 1990), chemistrytransport modeling (Vuolo et al., 2009), astrophysics (Almgren et al., 2010) etc.
The Després and Lagoutière (1999) scheme
The scheme of Després and Lagoutière (1999) is defined by their Eqs. (2) to (4). If ${F}_{i+\frac{\mathrm{1}}{\mathrm{2}}}>\mathrm{0}$, these equations brought to the notations of Eq. (12), give
with the same notations as for the Van Leer (1977) scheme (above). As above, if ($({\mathit{\alpha}}_{s,i}{\mathit{\alpha}}_{s,i\mathrm{1}})({\mathit{\alpha}}_{s,i+\mathrm{1}}{\mathit{\alpha}}_{s,i})\le \mathrm{0}$, no interpolation is performed and the scheme falls back to the simple Godunov donorcell formulation (Eq. 14). As stated by its authors, this scheme is antidiffusive. Unlike other schemes such as the Van Leer (1977) scheme described above, two unusual choices have been made by the authors in order to minimize diffusion by the advection scheme:

Their scheme is accurate only to the first order.

The scheme is linearly unstable, but nonlinearly stable (their Theorem 1).
The idea of the authors has been to make the interpolated value ${\stackrel{\mathrm{\u203e}}{\mathit{\alpha}}}_{s,i+\frac{\mathrm{1}}{\mathrm{2}}}$ as close as possible to the downstream value (${\mathit{\alpha}}_{s,i+\mathrm{1}}$ if the flux is positively oriented). This property is desirable because it is the key property in order to reduce numerical diffusion as much as mathematically possible while still maintaining the scheme stability. The authors present 1D case studies with their scheme obtaining extremely interesting results: fields that are initially concentrated on one single cell do not occupy more than three cells even after a long advection time (their Fig. 2), sharp gradients are very well preserved (their Fig. 1), and, more unexpectedly due to its antidiffusive character, the scheme also performs well in maintaining the shape of concentration fields with an initially smooth concentration gradient. After extensive testing, these authors also suggest (their Conjecture 1) that convergence of the simulated values towards exact values occur even if the time step is reduced before the space step: in simpler terms, this means that the scheme performs very well even at small CFL values, a property that is not shared by most advection schemes. Comparison of Eqs. (17) with (16) shows that the numerical cost of the Després and Lagoutière (1999) scheme is about the same as the Van Leer (1977) scheme.
2.1.2 2D discretization of advection and directional splitting
All the simulations in the present study are 2D x–z cases on a domain discretized over a regular Cartesian mesh. In order to work with the usual units (concentrations per unit volume and not per unit area), a third, degenerate space dimension has been introduced, with one grid cell in the y direction, δy=δx. Zero mass flux is prescribed in the y direction. Index i is attributed to the x direction, index k to the z direction, and no index is attributed to the degenerate y dimension. Eq. (12) becomes
Here, ${\stackrel{\mathrm{\u203e}}{F}}_{i\frac{\mathrm{1}}{\mathrm{2}},k}$ is the timeaveraged mass flux of air through the left boundary of cell i,k between t and t+Δt, and ${\stackrel{\mathrm{\u203e}}{F}}_{i+\frac{\mathrm{1}}{\mathrm{2}},k}$ and ${\stackrel{\mathrm{\u203e}}{F}}_{i,k\pm \frac{\mathrm{1}}{\mathrm{2}}}$ have similar definitions. ${\stackrel{\mathrm{\u203e}}{\mathit{\alpha}}}_{s,i\frac{\mathrm{1}}{\mathrm{2}},k}$ is the mixing ratio of species s in the air volume entering cell i,k through its left boundary between t and t+Δt. If 𝒱 is the geometric volume containing at time t all the air parcels that are going to cross the left boundary of cell i,k between t and t+Δt then
From a practical point of view, it is extremely difficult to actually find the contours of 𝒱 and to reconstruct and integrate the concentrations of air and of tracer over this volume. This is why, in practice, Eulerian models in Cartesian grids tend to split between the two (or three) space directions. Integrating separately direction x and then direction z over time Δt (generally called “Lie splitting”) gives order1 error, and applying the socalled Strang splitting (Strang, 1968) by integrating first in the x direction over $\frac{\mathrm{\Delta}t}{\mathrm{2}}$, then in the z direction over Δt, and finally once again in the x direction over $\frac{\mathrm{\Delta}t}{\mathrm{2}}$ reduces the splitting error to order 2.
2.1.3 Tested configurations
Since the present study is aimed at studying vertical transport only, we chose to test the Godunov, Van Leer (1977), Colella and Woodward (1984), and Després and Lagoutière (1999) schemes in the vertical direction with the same transport strategy over the x axis, namely the Colella and Woodward (1984) PPM scheme. For all simulations, splitting between the x and z direction has been performed using Strang splitting as described above, in order to maintain secondorder accuracy for the Van Leer and PPM simulations. While it would have been possible to use Lie splitting for simulations Upwind and DL99 (for which the vertical advection scheme is only firstorder accurate), the choice of using Strang splitting for all four simulations guarantees that all simulations are strictly identical except for their vertical advection scheme. The configurations that have been tested are summarized in Table 1.
Colella and Woodward (1984)Godunov and Bohachevsky (1959)Colella and Woodward (1984)Van Leer (1977)Colella and Woodward (1984)Colella and Woodward (1984)Colella and Woodward (1984)Després and Lagoutière (1999)2.2 Test case definition
We have defined three test cases designed to be representative of longrange tracer transport situations in the atmosphere. The simulation domain covers an x–z domain with length L=2000 km from west to east and thickness H=12 km, with periodic boundary conditions at the lateral boundaries and open boundaries at the top and at the bottom of the simulation domain, with clean air (zero tracer concentration) entering from these boundaries. We will use T=86 400 s, the length of a complete day on Earth, as the timescale for the case studies, along with the corresponding pulsation $\mathit{\omega}=\frac{\mathrm{2}\mathit{\pi}}{T}$. The number density of the carrying fluid (air) will be assumed uniform. These simplifying assumptions are designed to ease the formalism and the formulation of exact solutions of the problem. Two situations of relevance for atmospheric tracer transport have been represented, along with another numerical experiment designed in order to investigate the properties of the tested transport configurations in terms of convergence rate. Case 1, presented in Sect. 2.2.1, aims to represent the formation of a thin plume from an initially thicker tracer volume through the action of zonal wind shear, a situation typical of longrange advection of polluted plumes in the free troposphere. Case 2, presented in Sect. 2.2.2, represents longrange advection of a thin plume under the action of a strong zonal wind.
Except for tests of convergence rates for which increasingly fine discretizations have to be tested, the domain is discretized into 80 evenly spaced cells from west to east (Δx=25 km) and 24 evenly spaced cells from bottom to top (Δz=500 m).
2.2.1 Case 1: thin layer formation under wind shear
In this case, we consider the evolution of an inert tracer initially distributed as follows:
with h_{1}=1500 m the halfthickness of the initial layer, δx_{1}=25 km the halflength of the initial layer, H the height of domain top (H=12 km; see Table 2) and L the length of domain (L=2000 km). This describes a uniform plume initially confined vertically between z=4500 m and z=7500 m and horizontally between x=975 km and x=1025 km, with an initial (and arbitrary) mixing ratio of 100 ppb.
Zonal wind is constant in time, zonally uniform and vertically sheared:
with ${U}_{\mathrm{0}}=\frac{L}{\mathrm{2}T}$ so that the horizontal motion of fluid at $z=\frac{H}{\mathrm{2}}$ brings it back at its initial position after a time 2T, which will be the duration of the numerical experiment.
We add a vertical wind defined as follows:
The vertical wind speed scale is taken as $\mathrm{5}\times {\mathrm{10}}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}\mathrm{m}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{\mathrm{1}}$, a typical scale for synopticscale vertical motion in the troposphere. Since $\frac{\partial u}{\partial x}=\frac{\partial w}{\partial z}=\mathrm{0}$ and since the density field is uniform, this mass flux is nondivergent.
This case can describe a plume that is initially vertical, covering a 50 km wide column (two grid cells), uniformly spanned over a 3 km altitude range (4500 to 7500 m, corresponding to six grid cells; see Table 2). This initially thick vertical column later evolves under the action of wind shear into a thin layer.
Direct integration of Eq. (23) and then of Eq. (22) give access immediately to the position of a particle initially (t_{i}=0) at position (x_{i};z_{i}):
and
Equations (24) and (25) describe the superposition of a horizontal motion forced by the horizontal wind speed defined in Eq. (22), and an elliptic motion with pulsation ω due to the oscillating vertical speed (Eq. 23) and its interaction with the horizontal wind shear. The vertical semiax of this ellipse is $\frac{{w}_{\mathrm{0}}}{\mathit{\omega}}\simeq \mathrm{688}\phantom{\rule{0.125em}{0ex}}\mathrm{m}$, and the horizontal semiax is $\frac{\mathrm{2}{U}_{\mathrm{0}}{W}_{\mathrm{0}}}{H{\mathit{\omega}}^{\mathrm{2}}}\simeq \mathrm{9.12}\times {\mathrm{10}}^{\mathrm{3}}\phantom{\rule{0.125em}{0ex}}\mathrm{m}$.
With (x(t);z(t)) being, for any given time t, affine functions of (x_{i};z_{i}), the wind field of Eqs. (22) and (23) advects straight lines into straight lines. In particular, the initial rectangular zone containing the tracer will be advected, at any given time, into a parallelogram, whose summits are readily given by applying Eqs. (24) and (25) to the summits of the initial rectangle. Inside this moving parallelogram, that is increasingly tilted with time, tracer mixing ratio is equal to 100 ppb, zero outside, giving access to the exact solution of the case at any time.
2.2.2 Case 2: longrange advection of thin layer
In this case, the initial tracer mixing ratio is as follows:
with h_{2}=500 m. These equations describe a zonally infinite layer contained between z=5500 m and z=6500 m. As in Case 1, α_{m}=100 ppb.
Zonal wind is constant and uniform:
so that the horizontal motion of the fluid brings it back at its initial x coordinate after time 2T, which will be the duration of the experiment. Vertical wind is defined as follows:
As for Case 1, ${w}_{\mathrm{0}}=\mathrm{5}\times {\mathrm{10}}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}\mathrm{m}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{\mathrm{1}}$. This vertical wind speed has two maxima and two minima over the horizontal domain.
Integration of Eqs. (28) and (29) is immediate and give the trajectory of a particle located at (x_{i};z_{i}) at time t=0:
and
In particular, after a time kT, k being an integer, all the fluid particles are displaced by distance $\frac{kL}{\mathrm{2}}$ in the horizontal direction and back to their initial altitude. At these times, since the initial tracer plume is zonally uniform and infinite, the field of mixing ration will be exactly back to its initial value everywhere.
This case is a simplified representation of longrange advection of a 1 km thick layer of inert tracer under the action of a uniform zonal wind and variable vertical wind, representing for example in an extremely simplified way the advection of this layer through synopticscale structures. In the atmosphere, such fine layers of tracers are frequently formed by stretching of initially thicker polluted layers, as represented in Case 1.
2.2.3 Case 3: fine layer advection and convergence rate test
This case has been designed to study the numerical convergence rate of the various configurations that will be tested as a function of space resolution. The case setup is the same as for Case 2, with wind speeds similar to Eqs. (28) and (29):
and
Due to the need for increasing resolution and therefore the numerical cost of simulations, we simulate only one spatial and temporal period of this case (instead of two spatial and temporal periods for Case 2), hence the differences between Eqs. (32) and (33) and (28) and (29).
The initial tracer mixing ratio is prescribed as follows:
with h_{3}=1500 m. This initial mixing ratio distribution has been designed to be 𝒞^{2} (and in fact 𝒞^{3}) everywhere. It is therefore smooth enough to permit a convergence experiment with all the transport schemes which rely, at most, on the existence and continuity of the secondorder derivative of the transported field (for the PPM scheme).
For convergence rate tests, five different resolutions have been tested (Table 2).
2.3 Implementation
Implementation of the idealized experiments and transport scheme have been done within the under development ToyCTM code. ToyCTM is a Python code targeting chemistrytransport studies in academic cases. So far, ToyCTM relies on classical numpy arrays. Its objectoriented design provides a class structure enabling extensibility; i.e., users can easily code new transport schemes or define personal grid geometry. A basic chemistry module is present and allows one to test chemical reactions on top of transport. See the Code availability section at the end of the paper for access to the code version used for the present study and to the current development version of the code.
3.1 Case 1: thin layer formation under wind shear
Figure 1 shows the outputs of the four simulations realized with different schemes for vertical advection. All four simulations succeed in reproducing the tilted orientation of the final plume and its location but differ greatly in their maximal value and spatial extension, with the smallest maximal value of mixing ratio for the upwind scheme and the maximal value for the DL99 scheme, with the Van Leer and PPM schemes ranging in between. From a quantitative point of view (Table 3), the Upwind, VL and PPM schemes perform as could be expected from their order of accuracy, with the thirdorder PPM scheme performing better than the secondorder Van Leer scheme and the firstorder Upwind scheme in terms of all the diagnostics that have been calculated. More surprisingly, the firstorder DL99 scheme performs better than all these schemes in terms of all these diagnostics, by a wide margin: in this case study, the performance gain of DL99 relative to PPM is similar to the gain of PPM relative to Godunov in terms of maximal mixing ratio, percentage of mass in the envelope and accuracy.
3.2 Case 2: longrange advection of thin layer
Figure 2 shows the final state of Case 2 simulation for the four advection schemes that have been tested, as compared to the exact solution. Visually, the Després and Lagoutière (1999) scheme has performed best in bringing virtually all the tracer back into its original envelope after two complete vertical oscillations. Its performance in this case is the best of all the tested schemes (Table 4), with only a 7.4 % reduction in the maximal value of tracer mixing ratio (49 % with the PPM method, even more with the Godunov and Van Leer schemes), and very small error values in $\Vert \cdot {\Vert}_{\mathrm{1}}$ and $\Vert \cdot {\Vert}_{\mathrm{2}}$ compared to the other tested schemes. 90.6 % of the mass is contained in the theoretical envelope where it should be after the end of the numerical experiment (50.3 % only with the PPM scheme).
3.3 Case 3: fine layer advection and convergence rate test
The numerical convergence rates of all the tested advection configurations for $\Vert \cdot {\Vert}_{\mathrm{1}}$ and $\Vert \cdot {\Vert}_{\mathrm{2}}$ based on the last segment in Fig. 3a–b (between n_{x}=160 and n_{x}=320) are shown in Table 5. In $\Vert \cdot {\Vert}_{\mathrm{1}}$. These orders of convergence are around 0.8 for the DL99 and Upwind schemes because vertical resolution is still insufficient even at this high resolution to ensure theoretical convergence for these schemes. The Van Leer scheme yields a convergence rate of 1.80 in $\Vert \cdot {\Vert}_{\mathrm{1}}$, 2.43 for the PPM scheme. A smoother test case has been performed to check that all schemes are able to obtain convergence up to their theoretical order (see Appendix).
Figure 3 shows that, when resolution is too coarse to appropriately resolve the thin layer, the Després and Lagoutière (1999) scheme strongly overperforms higherorder schemes and offers an accuracy that is substantially better than both the Van Leer (1977) scheme and the Colella and Woodward (1984) scheme. This is consistent with the results obtained on Cases 1 and 2 and explains why in these cases the Després and Lagoutière (1999) scheme permits to obtain excellent results in reproducing the thin layer of tracer. On the other hand, when vertical resolution becomes sufficient to appropriately resolve the smooth tracer layer, higherorder schemes perform better than the Després and Lagoutière (1999) schemes, which in turns fall back towards an accuracy similar to the Godunov and Bohachevsky (1959) scheme, consistent with its expected order or accuracy.
When vertical resolution becomes much too fine compared to the size of the modeled object (polluted plume thicker that ≃50 Δz), accuracy of the simulations with the Després and Lagoutière (1999) scheme stops improving with resolution at an order1 rate (Fig. A1). Examination of the simulation outputs for these configurations reveal that progress in accuracy is hindered by undesirable smallscale oscillations that degrade accuracy (not shown).
The numerical experiments that have been presented confirm the interest of using the Després and Lagoutière (1999) scheme for vertical transport in chemistrytransport models, as already claimed by Lachatre et al. (2020). The idealized framework set up here permits examination of some of the questions that were out of reach in the real case of Lachatre et al. (2020) due to uncertainties in the forcing meteorological fields, the volcanic emissions and the lack of accurate measurements of plume structure.
The numerical experiments exposed here confirm that the Després and Lagoutière (1999) scheme is less diffusive than Van Leer (1977) and Colella and Woodward (1984), as could be expected due to its antidiffusive design. They also reveal that, in the presence of sharp vertical gradients that are not adequately resolved at model resolution, using the Després and Lagoutière (1999) scheme for vertical transport also increases model accuracy compared to the exact solution. This improvement is substantial for both Case 1 (Table 3 and Fig. 1) representing the formation of a thin polluted layer under the action of wind shear and Case 2 (Table 4 and Fig. 2) representing longrange advection of a thin polluted layer. The objective scores as well as the visual comparison of the simulated final state with the exact solution show that, at this resolution and for both these cases, using the Després and Lagoutière (1999) scheme reduces diffusion and increases accuracy compared to the schemes of Godunov and Bohachevsky (1959), Van Leer (1977), and Colella and Woodward (1984). While reduction in diffusion is in line with the results of Lachatre et al. (2020) and could be expected because of the design of the Després and Lagoutière (1999) schemes, improved accuracy in presence of sharp gradients is a strong argument in favor of using the Després and Lagoutière (1999) schemes for vertical transport in chemistrytransport models.
Improved accuracy of a loworder scheme compared to higherorder schemes for a given resolution is not impossible from a theoretical point of view but still counterintuitive since higherorder schemes are designed to reduce numerical error at any given resolution compared to lowerorder schemes due to “smarter” reconstruction procedures. Theory imposes that, if model resolution is fine enough and if the tracer field is smooth, higherorder schemes should be more accurate than lowerorder schemes. However, as shown by Godunov and Bohachevsky (1959), linear higherorder schemes cannot be monotonous, a property usually known as Godunov's theorem. This is why, to ensure monotonicity, the schemes of Van Leer (1977) and Colella and Woodward (1984) include nonlinear slopelimiters which are activated in the vicinity of extrema and discontinuities. In the vicinity of discontinuities, these formulations introduce large inaccuracies: in these schemes, the use of slopelimiters introduce large errors in the vicinity of discontinuities, and these errors generate excessive numerical diffusion, which is visible in Figs. 1 and 2. On the other hand, as discussed by its creators, the Després and Lagoutière (1999) scheme is designed to reduce numerical diffusion in these areas of steep gradients, which explains why it performs better than Van Leer (1977) and Colella and Woodward (1984) in all respects for Cases 1 and 2, which describe discontinuous tracer layers (Tables 3 and 4).
To understand this surprisingly good behavior of Després and Lagoutière (1999) compared to the higherorder Van Leer (1977) and Colella and Woodward (1984) schemes, we have performed a convergence test for advection of a 3000 m thick layer with a smooth (𝒞^{3}) initial profile for the tracer mixing ratio (Eq. 34). This convergence test (Fig. 3) shows that the Després and Lagoutière (1999) scheme performs better than these classical order2 schemes if model vertical resolution Δz is equal to 1000 or 500 m, but that due to their faster convergence rate, order2 schemes perform better if Δz≤250 m. At these coarse resolutions, gradients in the initial tracer field (Eq. 34) appear so sharp that the Després and Lagoutière (1999) scheme yields improved accuracy compared to Van Leer (1977) and Colella and Woodward (1984) due to its better ability to deal with sharp gradients without introducing excessive numerical diffusion. On the other hand, when resolution is fine enough, the tracer field and its successive derivatives do not vary drastically from one model cell to the next, so that the linear interpolations of Van Leer (1977) and the quadratic interpolations of Colella and Woodward (1984) work adequately and reduce numerical errors compared to the firstorder evaluation of concentration of Després and Lagoutière (1999).
In more general words, this result suggests that the Després and Lagoutière (1999) scheme may be expected to perform better than classical schemes in chemistrytransport models for the advection of polluted plumes thinner than ≃6 Δz (Δz being the model's vertical resolution), while higherorder schemes can be expected to perform better for the advection of polluted plumes thicker than 6 Δz if we suppose that the plume has a smooth vertical profile. Under realistic conditions of wind shear, these conditions of sufficient smoothness and thickness might actually be very difficult to reach since, as described in Case 1, vertical wind shear tend to the permanent thinning of atmospheric plumes (this question is discussed in detail in Zhuang et al., 2018) so that the Després and Lagoutière (1999) may frequently overperform classical order2 schemes in realistic wind conditions including wind shear.
The firstorder, antidiffusive advection scheme of Després and Lagoutière (1999) has been compared to classical secondorder schemes of Van Leer (1977) and Colella and Woodward (1984) in three idealized test cases representing longrange atmospheric transport of thin plumes in the troposphere. It is shown that the Després and Lagoutière (1999) generally overperforms these two schemes, offering both improved accuracy and reduced diffusion. The Després and Lagoutière (1999) scheme is shown to allow a correct representation of longrange advection of fine tracer layers when vertical resolution is so coarse that the classical Van Leer and PPM schemes exhibit much weaker performance due to excessive vertical diffusion, a feature that appears in the present idealized test cases but also in realistic cases (e.g., Colette et al., 2010; Lachatre et al., 2020).
Convergence tests show that this improved performance exists only for tracer layers that are represented by less than six grid cells in the vertical direction, while for finer resolutions and smooth initial tracer fields, higherorder schemes perform better than the Després and Lagoutière (1999) firstorder scheme due to their faster convergence towards the exact solution. This suggests that, if model resolution is fine enough to represent properly the tracer plumes and their concentration gradients, higherorder schemes may still be a better choice.
We think that these results are important because they explain under which conditions the Després and Lagoutière (1999) is able to reduce excessive vertical diffusion in CTMs, as was observed in Lachatre et al. (2020), and that this reduced vertical diffusion comes together with an improved accuracy compared to the exact solution. Since the vertical resolution of chemistrytransport models is usually coarse in the free troposphere, while the advected plumes tend to be extremely thin in this part of the atmosphere due to the action of wind shear, the present study along with Lachatre et al. (2020) advocates for using the Després and Lagoutière (1999) transport scheme for chemistrytransport modeling in the free troposphere, and probably even more in the stratosphere where vertical diffusion needs to be extremely reduced. It is also worth noting that Després and Lagoutière (1999) have shown that their scheme maintains its convergence and lowdiffusion properties even if the CFL number becomes small, which is very common for vertical advection in the free troposphere due to the typically small vertical speed of air motion (typically a few centimeters per second).
More investigation is needed in real and/or idealized cases to address several questions:

Does the Després and Lagoutière (1999) perform better in the boundary layer, where CTM resolution is typically finer than in the free troposphere?

If not, is it possible to use a traditional transport scheme in the boundary layer and the Després and Lagoutière (1999) scheme in the free troposphere without introducing numerical artifacts in the buffer zone?

Atmospheric chemistry being a nonlinear process, how does reduction in excessive numerical diffusion in the troposphere affect representation of chemistry inside chemically active air masses such as volcanic or biomass burning plumes?

Is it desirable to use antidiffusive transport schemes in the horizontal directions as well, and under which conditions should they be used?
A1 Case definition
This case has been designed to study the numerical convergence rate of the various configurations that will be tested as a function of space resolution. The case setup is the same as for Cases 2 and 3, with wind speeds from Eqs. (28) and (29), but the initial tracer mixing ratio is prescribed as follows:
with ${h}_{\mathrm{4}}=\frac{H}{\mathrm{2}}=\mathrm{6000}\phantom{\rule{0.125em}{0ex}}\mathrm{m}$ and $\mathit{\delta}=\frac{L}{\mathrm{2}}=\mathrm{5}\times {\mathrm{10}}^{\mathrm{5}}\phantom{\rule{0.125em}{0ex}}\mathrm{m}$. This initial tracer distribution represents a cosine bell centered at $(x=\frac{L}{\mathrm{2}};z=\frac{H}{\mathrm{2}})$, C^{∞} everywhere with extremely smooth variations in space. Domain resolution and sizes are as shown in Table 2.
A2 Results
Table A1 and Fig. A1 show that all the tested configurations exhibit the expected convergence order at least in $\Vert \cdot {\Vert}_{\mathrm{1}}$, but accuracy with the DL99 scheme stops improving when vertical resolution becomes “too fine” compared to the thickness of the represented layer. In the present case, this “saturation” of convergence occurs between nx=80 (nz=48) and nx=160 (nz=96):
ToyCTM is free software distributed under the GNU General Public License v2. The exact version used for the present study (Mailler and Pennel, 2020) is available permanently in the HAL repository at https://hal.archivesouvertes.fr/hal02933095. The latest stable version of ToyCTM is available from https://gitlab.in2p3.fr/ipsl/lmd/intro/toyctm//archive/master/toyctmmaster.tar.gz (last access: 29 April 2021).
All the authors have contributed to the design of the simulated cases. SM performed and analyzed the simulations and developed the software with RP. LM, ML and RP contributed to writing and improving the paper.
The authors declare that they have no conflict of interest.
The simulations have been performed at the ESPRI/IPSL data center and at TGCC under GENCI A0070110274 allocation.
This research has been supported by the Agence de l'Innovation de Défense (TROMPET grant).
This paper was edited by Simone Marras and reviewed by two anonymous referees.
Almgren, A. S., Beckner, V. E., Bell, J. B., Day, M. S., Howell, L. H., Joggerst, C. C., Lijewski, M. J., Nonaka, A., Singer, M., and Zingale, M.: CASTRO: A New Compressible Astrophysical Solver. I. Hydrodynamics and Selfgravity, Astrophys. J., 715, 1221–1238, https://doi.org/10.1088/0004637X/715/2/1221, 2010. a, b
Carpenter Jr., R. L., Droegemeier, K. K., Woodward, P. R., and Hane, C. E.: Application of the Piecewise Parabolic Method (PPM) to Meteorological Modeling, Mon. Weather Rev., 118, 586–612, https://doi.org/10.1175/15200493(1990)118<0586:AOTPPM>2.0.CO;2, 1990. a
Colella, P. and Sekora, M. D.: A limiter for PPM that preserves accuracy at smooth extrema, J. Comput. Phys., 227, 7069–7076, https://doi.org/10.1016/j.jcp.2008.03.034, 2008. a
Colella, P. and Woodward, P. R.: The piecewise parabolic method (PPM) for gasdynamical simulations, J. Comput. Phys., 11, 38–39, 1984. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x
Colette, A., Alsac, N., Bessagnet, B., Biaudet, H., Chiappini, L., Favez, O., Frejafon, E., Gautier, F., Godefroy, F., Haeffelin, M., Leoz, E., Malherbe, L., Meleux, F., Menut, L., Morille, Y., Papin, A., Pietras, C., Ramel, M., and Rouil, L.: Assessment of the impact of the Eyjafjallajökull's eruption on surface air quality in France, Atmos. Environ., 12, 1217–1221, https://doi.org/10.1016/j.atmosenv.2010.09.064, 2010. a, b, c
Després, B. and Lagoutière, F.: Un schéma non linéaire antidissipatif pour l'équation d'advection linéaire, CR l'Acad. Sci. IMath., 328, 939–943, https://doi.org/10.1016/S07644442(99)803012, 1999. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x, y, z, aa, ab, ac, ad, ae, af, ag, ah, ai, aj, ak, al
Eastham, S. D. and Jacob, D. J.: Limits on the ability of global Eulerian models to resolve intercontinental transport of chemical plumes, Atmos. Chem. Phys., 17, 2543–2553, https://doi.org/10.5194/acp1725432017, 2017. a
Emery, C., Tai, E., Yarwood, G., and Morris, R.: Investigation into approaches to reduce excessive vertical transport over complex terrain in a regional photochemical grid model, Atmos. Environ., 45, 7341–7351, https://doi.org/10.1016/j.atmosenv.2011.07.052, 2011. a, b, c, d
Godunov, S. K. and Bohachevsky, I.: Finite difference method for numerical computation of discontinuous solutions of the equations of fluid dynamics, Matematičeskij sbornik, 47, 271–306, available at: https://hal.archivesouvertes.fr/hal01620642 (last access: 27 April 2021), 1959. a, b, c, d, e, f
Jöckel, P., von Kuhlmann, R., Lawrence, M. G., Steil, B., Brenninkmeijer, C. A. M., Crutzen, P. J., Rasch, P. J., and Eaton, B.: On a fundamental problem in implementing fluxform advection schemes for tracer transport in 3dimensional general circulation and chemistry transport models, Q. J. Roy. Meteor. Soc., 127, 1035–1052, https://doi.org/10.1002/qj.49712757318, 2001. a
Lachatre, M., Mailler, S., Menut, L., Turquety, S., Sellitto, P., Guermazi, H., Salerno, G., Caltabiano, T., and Carboni, E.: New strategies for vertical transport in chemistry transport models: application to the case of the Mount Etna eruption on 18 March 2012 with CHIMERE v2017r4, Geosci. Model Dev., 13, 5707–5723, https://doi.org/10.5194/gmd1357072020, 2020. a, b, c, d, e, f, g, h, i, j, k
Mailler, S. and Pennel, R.: toyCTM, available at: https://hal.archivesouvertes.fr/hal02933095 (last access: 27 April 2021), 2020. a
Mailler, S., Menut, L., Khvorostyanov, D., Valari, M., Couvidat, F., Siour, G., Turquety, S., Briant, R., Tuccella, P., Bessagnet, B., Colette, A., Létinois, L., Markakis, K., and Meleux, F.: CHIMERE2017: from urban to hemispheric chemistrytransport modeling, Geosci. Model Dev., 10, 2397–2423, https://doi.org/10.5194/gmd1023972017, 2017. a, b, c, d
Strang, G.: On the construction and comparison of difference schemes, SIAM J. Numer. Anal., 5, 506–517, https://doi.org/10.1137/0705041, 1968. a
Van Leer, B.: Towards the ultimate conservative difference scheme. IV. A new approach to numerical convection, J. Comput. Phys., 23, 276–299, https://doi.org/10.1016/00219991(77)90095X, 1977. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r
Vuolo, M. R., Menut, L., and Chepfer, H.: Impact of transport schemes on modeled dust concentrations, J. Atmos. Ocean. Tech., 26, 1135–1143, 2009. a, b
Zhuang, J., Jacob, D. J., and Eastham, S. D.: The importance of vertical resolution in the free troposphere for modeling intercontinental plumes, Atmos. Chem. Phys., 18, 6039–6055, https://doi.org/10.5194/acp1860392018, 2018. a, b