Using the Després and Lagoutière (1999) antidiffusive transport scheme: a promising and novel method against excessive vertical diffusion in chemistry-transport models

The potential of the antidiffusive transport scheme proposed by Després and Lagoutière (1999) for resolving vertical transport in chemistry-transport models is investigated in an idealized framework with very encouraging results. We show that, compared to classical higher-order 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 chemistry-transport 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 chemistry-transport models. Therefore, we think that this scheme is an extremely promising solution for reducing numerical diffusion in chemistry-transport models.

1 Introduction

Reducing numerical errors in chemistry-transport models (CTMs) is a necessary task for those wishing to improve these models' performance. Among the well-known errors in Eulerian CTMs, excessive numerical diffusion in all directions is a well-known 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 ground-level 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 fine-scale 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 model-specific 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 first-order to a second-order 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 long-range 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 chemistry-transport 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 long-range 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 chemistry-transport model fed by real-life 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 above-cited schemes and to the order-1 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.

2 Numerical methods and case description

Continuity equation for the motion of air is as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle \frac{\partial C}{\partial t}}+\mathrm{\nabla }\left(\mathbf{\Phi }\right)=\mathrm{0},\end{array}$$

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:

$$\begin{array}{}\text{(2)}& {\displaystyle }{\displaystyle \frac{\partial C_s}{\partial t}}+\mathrm{\nabla }\left({\mathbf{\Phi }}_s\right)=\mathrm{0},\end{array}$$

where C_s is concentration of species s in molecules per unit volume and Φ_s=C_su is the flux vector for species s. Equivalently, Eq. (2) becomes

$$\begin{array}{}\text{(3)}& {\displaystyle }{\displaystyle \frac{\partial C_s}{\partial t}}+\mathrm{\nabla }\left({\mathit{\alpha }}_s\mathbf{\Phi }\right)=\mathrm{0},\end{array}$$

where α_s is the mixing ratio for species s:

$$\begin{array}{}\text{(4)}& {\displaystyle }{\displaystyle }{\mathit{\alpha }}_s={\displaystyle \frac{C_s}{C}}.\end{array}$$

Chemistry-transport 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

$$\begin{array}{}\text{(5)}& {\displaystyle }{\displaystyle \frac{\partial C_s}{\partial t}}+{\displaystyle \frac{\partial \left(C_s\mathit{u}\right)}{\partial x}}=\mathrm{0}.\end{array}$$

Here we follow a semi-Lagrangian 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{‾}}{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{‾}}{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:

$$\begin{array}{}\text{(6)}& {\displaystyle }{\displaystyle \frac{C_{s,i}\left(t+\mathrm{\Delta }t\right)-C_{s,i}\left(t\right)}{\mathrm{\Delta }t}}={\displaystyle \frac{{\stackrel{\mathrm{‾}}{u}}_{i-\frac{\mathrm{1}}{\mathrm{2}}}{\stackrel{\mathrm{‾}}{C}}_{s,i-\frac{\mathrm{1}}{\mathrm{2}}}-{\stackrel{\mathrm{‾}}{u}}_{i+\frac{\mathrm{1}}{\mathrm{2}}}{\stackrel{\mathrm{‾}}{C}}_{s,i+\frac{\mathrm{1}}{\mathrm{2}}}}{x_{i+\frac{\mathrm{1}}{\mathrm{2}}}-x_{i-\frac{\mathrm{1}}{\mathrm{2}}}}},\end{array}$$

where

$$\begin{array}{}\text{(7)}& {\displaystyle }{\displaystyle }{\stackrel{\mathrm{‾}}{C}}_{s,i-\frac{\mathrm{1}}{\mathrm{2}}}={\displaystyle \frac{\mathrm{1}}{{\mathrm{\Delta }}^{-}x}}\underset{x_{i-\frac{\mathrm{1}}{\mathrm{2}}}-{\mathrm{\Delta }}^{-}x}{\overset{x_{i-\frac{\mathrm{1}}{\mathrm{2}}}}{\int }}{C}_s(x,t)\mathrm{d}x\end{array}$$

and

$$\begin{array}{}\text{(8)}& {\displaystyle }{\displaystyle }{\stackrel{\mathrm{‾}}{C}}_{s,i+\frac{\mathrm{1}}{\mathrm{2}}}={\displaystyle \frac{\mathrm{1}}{{\mathrm{\Delta }}^{+}x}}\underset{x_{i+\frac{\mathrm{1}}{\mathrm{2}}}}{\overset{x_{i+\frac{\mathrm{1}}{\mathrm{2}}}+{\mathrm{\Delta }}^{+}x}{\int }}{C}_s(x,t)\mathrm{d}x.\end{array}$$

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:

$$\begin{array}{}\text{(9)}& {\displaystyle }{\displaystyle \frac{C_i(t+\mathrm{\Delta }t)-C_i\left(t\right)}{\mathrm{\Delta }t}}={\displaystyle \frac{{\stackrel{\mathrm{‾}}{u}}_{i-\frac{\mathrm{1}}{\mathrm{2}}}{\stackrel{\mathrm{‾}}{C}}_{i-\frac{\mathrm{1}}{\mathrm{2}}}-{\stackrel{\mathrm{‾}}{u}}_{i+\frac{\mathrm{1}}{\mathrm{2}}}{\stackrel{\mathrm{‾}}{C}}_{i+\frac{\mathrm{1}}{\mathrm{2}}}}{x_{i+\frac{\mathrm{1}}{\mathrm{2}}}-x_{i-\frac{\mathrm{1}}{\mathrm{2}}}}}.\end{array}$$

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

$$\begin{array}{}\text{(10)}& {\displaystyle }{\displaystyle }{\stackrel{\mathrm{‾}}{\mathit{\alpha }}}_{s,i\pm \frac{\mathrm{1}}{\mathrm{2}}}={\displaystyle \frac{{\stackrel{\mathrm{‾}}{C}}_{s,i\pm \frac{\mathrm{1}}{\mathrm{2}}}}{{\stackrel{\mathrm{‾}}{C}}_{i\pm \frac{\mathrm{1}}{\mathrm{2}}}}}\end{array}$$

and

$$\begin{array}{}\text{(11)}& {\displaystyle }{\displaystyle }{\stackrel{\mathrm{‾}}{F}}_{i\pm \frac{\mathrm{1}}{\mathrm{2}}}={\stackrel{\mathrm{‾}}{C}}_{i\pm \frac{\mathrm{1}}{\mathrm{2}}}{\stackrel{\mathrm{‾}}{u}}_{s,i\pm \frac{\mathrm{1}}{\mathrm{2}}}.\end{array}$$

With these notations, Eqs. (6) and (9) become

$$\begin{array}{}\text{(12)}& {\displaystyle }{\displaystyle \frac{C_{s,i}(t+\mathrm{\Delta }t)-C_{s,i}\left(t\right)}{\mathrm{\Delta }t}}={\displaystyle \frac{{\stackrel{\mathrm{‾}}{F}}_{i-\frac{\mathrm{1}}{\mathrm{2}}}{\stackrel{\mathrm{‾}}{\mathit{\alpha }}}_{s,i-\frac{\mathrm{1}}{\mathrm{2}}}-{\stackrel{\mathrm{‾}}{F}}_{i+\frac{\mathrm{1}}{\mathrm{2}}}{\stackrel{\mathrm{‾}}{\mathit{\alpha }}}_{s,i+\frac{\mathrm{1}}{\mathrm{2}}}}{x_{i+\frac{\mathrm{1}}{\mathrm{2}}}-x_{i-\frac{\mathrm{1}}{\mathrm{2}}}}}\end{array}$$

and

$$\begin{array}{}\text{(13)}& {\displaystyle }{\displaystyle \frac{C_i(t+\mathrm{\Delta }t)-C_i\left(t\right)}{\mathrm{\Delta }t}}={\displaystyle \frac{{\stackrel{\mathrm{‾}}{F}}_{i-\frac{\mathrm{1}}{\mathrm{2}}}-{\stackrel{\mathrm{‾}}{F}}_{i+\frac{\mathrm{1}}{\mathrm{2}}}}{x_{i+\frac{\mathrm{1}}{\mathrm{2}}}-x_{i-\frac{\mathrm{1}}{\mathrm{2}}}}}.\end{array}$$

Equations (12) and (13) are a flux-form reformulation of semi-Lagrangian 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{‾}}{\mathit{\alpha }}}_{s,i-\frac{\mathrm{1}}{\mathrm{2}}}$ lies between ${\stackrel{\mathrm{‾}}{\mathit{\alpha }}}_{s,i-\mathrm{1}}$ and ${\stackrel{\mathrm{‾}}{\mathit{\alpha }}}_{s,i}$ (and if ${\stackrel{\mathrm{‾}}{\mathit{\alpha }}}_{s,i+\frac{\mathrm{1}}{\mathrm{2}}}$ lies between ${\stackrel{\mathrm{‾}}{\mathit{\alpha }}}_{s,i}$ and ${\stackrel{\mathrm{‾}}{\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 chemistry-transport models usually resolve advection of trace species using an approach based on Eq. (12) rather than a straightforward resolution of Eq. (6).

Regarding chemistry-transport models, in practice, approximated values of ${\stackrel{\mathrm{‾}}{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 mass-wind inconsistencies discussed in, for example, Jöckel et al. (2001), Emery et al. (2011) and Lachatre et al. (2020), by working with analytically defined non-divergent 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 chemistry-transport model whose task, in flux-form, is to estimate the values of ${\stackrel{\mathrm{‾}}{\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 donor-cell scheme

The most simple way of estimating ${\stackrel{\mathrm{‾}}{\mathit{\alpha }}}_{s,i\pm \frac{\mathrm{1}}{\mathrm{2}}}$ is the Godunov donor-cell scheme (adapted from Godunov and Bohachevsky, 1959), simply evaluating ${\stackrel{\mathrm{‾}}{\mathit{\alpha }}}_{s,i,k+\frac{\mathrm{1}}{\mathrm{2}}}$ as follows:

$$\begin{array}{}\text{(14)}& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{‾}}{\mathit{\alpha }}}_{s,i+\frac{\mathrm{1}}{\mathrm{2}}}={\mathit{\alpha }}_{s,i}\text{if }{\stackrel{\mathrm{‾}}{F}}_{i+\frac{\mathrm{1}}{\mathrm{2}}}\ge \mathrm{0},\text{(15)}& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{‾}}{\mathit{\alpha }}}_{s,i+\frac{\mathrm{1}}{\mathrm{2}}}={\mathit{\alpha }}_{s,i+\mathrm{1}}\text{if }{\stackrel{\mathrm{‾}}{F}}_{i+\frac{\mathrm{1}}{\mathrm{2}}}<\mathrm{0}.\end{array}$$

This order-1 scheme is cheap, robust, linear, monotonous and mass-conservative but extremely diffusive. It is therefore important to find more accurate ways to estimate ${\stackrel{\mathrm{‾}}{\mathit{\alpha }}}_{s,i+\frac{\mathrm{1}}{\mathrm{2}}}$.

The Van Leer (1977) scheme

The second-order slope-limited scheme of Van Leer (1977) brought to our notations and assuming uniform air concentration yields the following expression of ${\stackrel{\mathrm{‾}}{\mathit{\alpha }}}_{s,i+\frac{\mathrm{1}}{\mathrm{2}}}$ (for ${\stackrel{\mathrm{‾}}{F}}_{i+\frac{\mathrm{1}}{\mathrm{2}}}>\mathrm{0}$):

$$\begin{array}{ll}{\displaystyle }{\stackrel{\mathrm{‾}}{\mathit{\alpha }}}_{s,i+\frac{\mathrm{1}}{\mathrm{2}}}=& {\displaystyle }{\mathit{\alpha }}_{s,i}+{\displaystyle \frac{\mathrm{1}-\mathit{\nu }}{\mathrm{2}}}\text{sign}\left({\mathit{\alpha }}_{s,i+\mathrm{1}}-{\mathit{\alpha }}_{s,i}\right)\\ {\displaystyle }& {\displaystyle }\times \text{Min}({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left|{\mathit{\alpha }}_{s,i+\mathrm{1}}-{\mathit{\alpha }}_{s,i-\mathrm{1}}\right|,\mathrm{2}\left|{\mathit{\alpha }}_{s,i+\mathrm{1}}-{\mathit{\alpha }}_{s,i}\right|,\\ \text{(16)}& {\displaystyle }& {\displaystyle }\mathrm{2}\left|{\mathit{\alpha }}_{s,i}-{\mathit{\alpha }}_{s,i-\mathrm{1}}\right|),\end{array}$$

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{‾}}{\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 donor-cell formulation (Eq. 14). This order-2 scheme has been used for decades in chemistry-transport modeling, being a good tradeoff between reasonably weak diffusion, at least compared to more simple schemes such as the Godunov donor-cell scheme, and small computational burden compared to higher-order 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 third-order by design, application of limiters in the vicinity of extrema introduces first-order truncation errors in their vicinity so that third-order 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 order-2 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), chemistry-transport 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

$$\begin{array}{ll}{\displaystyle }{\stackrel{\mathrm{‾}}{\mathit{\alpha }}}_{s,i\frac{\mathrm{1}}{\mathrm{2}}}=& {\displaystyle }{\mathit{\alpha }}_{s,k}+{\displaystyle \frac{\mathrm{1}-\mathit{\nu }}{\mathrm{2}}}\\ {\displaystyle }& {\displaystyle }\times \text{Max}\left[\mathrm{0},\text{Min}\left({\displaystyle \frac{\mathrm{2}}{\mathit{\nu }}}{\displaystyle \frac{{\mathit{\alpha }}_{s,i}-{\mathit{\alpha }}_{s,i-\mathrm{1}}}{{\mathit{\alpha }}_{s,i+\mathrm{1}}-{\mathit{\alpha }}_{s,i}}},{\displaystyle \frac{\mathrm{2}}{\mathrm{1}-\mathit{\nu }}}\right)\right]\\ \text{(17)}& {\displaystyle }& {\displaystyle }\times \left({\mathit{\alpha }}_{s,i+\mathrm{1}}-{\mathit{\alpha }}_{s,i}\right),\end{array}$$

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 donor-cell 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 non-linearly stable (their Theorem 1).

The idea of the authors has been to make the interpolated value ${\stackrel{\mathrm{‾}}{\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

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{C_{s,i,k}\left(t+\mathrm{\Delta }t\right)-C_{s,i,k}\left(t\right)}{\mathrm{\Delta }t}}\\ \text{(18)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{{\stackrel{\mathrm{‾}}{F}}_{i-\frac{\mathrm{1}}{\mathrm{2}},k}{\stackrel{\mathrm{‾}}{\mathit{\alpha }}}_{s,i-\frac{\mathrm{1}}{\mathrm{2}},k}-{\stackrel{\mathrm{‾}}{F}}_{i+\frac{\mathrm{1}}{\mathrm{2}},k}{\stackrel{\mathrm{‾}}{\mathit{\alpha }}}_{s,i+\frac{\mathrm{1}}{\mathrm{2}},k}}{x_{i+\frac{\mathrm{1}}{\mathrm{2}}}-x_{i-\frac{\mathrm{1}}{\mathrm{2}}}}}{\displaystyle }& {\displaystyle }+{\displaystyle \frac{{\stackrel{\mathrm{‾}}{F}}_{i,k-\frac{\mathrm{1}}{\mathrm{2}}}{\stackrel{\mathrm{‾}}{\mathit{\alpha }}}_{s,i,k-\frac{\mathrm{1}}{\mathrm{2}}}-{\stackrel{\mathrm{‾}}{F}}_{i,k+\frac{\mathrm{1}}{\mathrm{2}}}{\stackrel{\mathrm{‾}}{\mathit{\alpha }}}_{s,i,k+\frac{\mathrm{1}}{\mathrm{2}}}}{z_{k+\frac{\mathrm{1}}{\mathrm{2}}}-z_{k-\frac{\mathrm{1}}{\mathrm{2}}}}}.\end{array}$$

Here, ${\stackrel{\mathrm{‾}}{F}}_{i-\frac{\mathrm{1}}{\mathrm{2}},k}$ is the time-averaged mass flux of air through the left boundary of cell i,k between t and t+Δt, and ${\stackrel{\mathrm{‾}}{F}}_{i+\frac{\mathrm{1}}{\mathrm{2}},k}$ and ${\stackrel{\mathrm{‾}}{F}}_{i,k\pm \frac{\mathrm{1}}{\mathrm{2}}}$ have similar definitions. ${\stackrel{\mathrm{‾}}{\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

$$\begin{array}{}\text{(19)}& {\displaystyle }{\displaystyle }{\stackrel{\mathrm{‾}}{\mathit{\alpha }}}_{s,i-\frac{\mathrm{1}}{\mathrm{2}},k}={\displaystyle \frac{{\int }_{\mathcal{V}}{C}_s\left(x,z,t\right)\mathrm{d}\mathcal{V}}{{\int }_{\mathcal{V}}C\left(x,z,t\right)\mathrm{d}\mathcal{V}}}.\end{array}$$

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 order-1 error, and applying the so-called 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 second-order 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 first-order 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)

Table 1 Summary of the different transport configurations that have been tested.

Table 2 Domain resolution and size for Cases 1 and 2; set of increasing resolutions used for convergence tests (Cases 3 and 4). For all configurations, $nz*\mathrm{\Delta }z=\mathrm{12}\mathrm{000}\mathrm{m}$ is the domain vertical extension.

2.2 Test case definition

We have defined three test cases designed to be representative of long-range 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 long-range advection of polluted plumes in the free troposphere. Case 2, presented in Sect. 2.2.2, represents long-range 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:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathit{\alpha }\left(t=\mathrm{0},x,z\right)={\mathit{\alpha }}_m\text{ if }{\displaystyle \frac{H}{\mathrm{2}}}-h_{\mathrm{1}}\le z\le {\displaystyle \frac{H}{\mathrm{2}}}+h_{\mathrm{1}}\text{ and }\\ \text{(20)}& {\displaystyle }& {\displaystyle }{\displaystyle \frac{L}{\mathrm{2}}}-\mathit{\delta }{x}_{\mathrm{1}}\le x\le {\displaystyle \frac{L}{\mathrm{2}}}+\mathit{\delta }{x}_{\mathrm{1}},\text{(21)}& {\displaystyle }& {\displaystyle }\mathit{\alpha }\left(t=\mathrm{0},x,z\right)=\mathrm{0}\text{ otherwise,}\end{array}$$

with h₁=1500 m the half-thickness of the initial layer, δx₁=25 km the half-length 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:

$$\begin{array}{}\text{(22)}& {\displaystyle }{\displaystyle }u(x,z,t)=U_{\mathrm{0}}{\displaystyle \frac{\mathrm{2}z}{H}},\end{array}$$

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:

$$\begin{array}{}\text{(23)}& {\displaystyle }{\displaystyle }w(x,y,z,t)=w_{\mathrm{0}}\text{cos}\left(\mathit{\omega }T\right).\end{array}$$

The vertical wind speed scale is taken as $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{2}}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$, a typical scale for synoptic-scale 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 non-divergent.

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):

$$\begin{array}{}\text{(24)}& {\displaystyle }{\displaystyle }x\left(t\right)=x_i+{\displaystyle \frac{\mathrm{2}{U}_{\mathrm{0}}}{H}}{z}_{i}t+{\displaystyle \frac{\mathrm{2}{U}_{\mathrm{0}}{W}_{\mathrm{0}}}{H{\mathit{\omega }}^{\mathrm{2}}}}\left(\mathrm{1}-\text{cos}\left(\mathit{\omega }t\right)\right)\end{array}$$

and

$$\begin{array}{}\text{(25)}& {\displaystyle }{\displaystyle }z\left(t\right)=z_i+{\displaystyle \frac{W_{\mathrm{0}}}{\mathit{\omega }}}\text{sin}\left(\mathit{\omega }t\right).\end{array}$$

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 semi-ax of this ellipse is $\frac{w_{\mathrm{0}}}{\mathit{\omega }}\simeq \mathrm{688}\mathrm{m}$, and the horizontal semi-ax is $\frac{\mathrm{2}{U}_{\mathrm{0}}{W}_{\mathrm{0}}}{H{\mathit{\omega }}^{\mathrm{2}}}\simeq \mathrm{9.12}\times {\mathrm{10}}^{\mathrm{3}}\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: long-range advection of thin layer

In this case, the initial tracer mixing ratio is as follows:

$$\begin{array}{}\text{(26)}& {\displaystyle }& {\displaystyle }\mathit{\alpha }(t=\mathrm{0},x,z)={\mathit{\alpha }}_m\text{ if }-h_{\mathrm{2}}\le z-{\displaystyle \frac{H}{\mathrm{2}}}\le h_{\mathrm{2}},\text{(27)}& {\displaystyle }& {\displaystyle }\mathit{\alpha }(t=\mathrm{0},x,z)=\mathrm{0}\text{ otherwise,}\end{array}$$

with h₂=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:

$$\begin{array}{}\text{(28)}& {\displaystyle }{\displaystyle }u(x,z,t)=U_{\mathrm{0}}={\displaystyle \frac{L}{\mathrm{2}T}}\end{array}$$

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:

$$\begin{array}{}\text{(29)}& {\displaystyle }{\displaystyle }w(x,y,z,t)=w_{\mathrm{0}}\text{cos}\left({\displaystyle \frac{\mathrm{4}\mathit{\pi }x}{L}}\right).\end{array}$$

As for Case 1, $w_{\mathrm{0}}=\mathrm{5}\times {\mathrm{10}}^{-\mathrm{2}}\mathrm{m}{\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:

$$\begin{array}{}\text{(30)}& {\displaystyle }{\displaystyle }x\left(t\right)={\displaystyle \frac{Lt}{\mathrm{2}T}}\end{array}$$

and

$$\begin{array}{}\text{(31)}& {\displaystyle }{\displaystyle }z\left(t\right)=z_i+{\displaystyle \frac{W_{\mathrm{0}}T}{\mathrm{2}\mathit{\pi }}}\left(\text{sin}\left({\displaystyle \frac{\mathrm{4}\mathit{\pi }{x}_i}{L}}+{\displaystyle \frac{\mathrm{2}\mathit{\pi }t}{T}}\right)-\text{sin}\left({\displaystyle \frac{\mathrm{4}\mathit{\pi }{x}_i}{L}}\right)\right).\end{array}$$

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 long-range 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 synoptic-scale 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):

$$\begin{array}{}\text{(32)}& {\displaystyle }{\displaystyle }u(x,z,t)=U_{\mathrm{0}}={\displaystyle \frac{L}{T}}\end{array}$$

and

$$\begin{array}{}\text{(33)}& {\displaystyle }{\displaystyle }w(x,y,z,t)=w_{\mathrm{0}}\text{cos}\left({\displaystyle \frac{\mathrm{2}\mathit{\pi }x}{L}}\right).\end{array}$$

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:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{\alpha }}_i(x,z)={\displaystyle \frac{{\mathit{\alpha }}_m}{\mathrm{4}}}{\left(\mathrm{1}+\text{cos}\left({\displaystyle \frac{\mathit{\pi }(z-H/\mathrm{2})}{h_{\mathrm{3}}}}\right)\right)}^{\mathrm{2}}\\ \text{(34)}& {\displaystyle }& {\displaystyle }\text{ if }-h_{\mathrm{3}}\le z-{\displaystyle \frac{H}{\mathrm{2}}}\le h_{\mathrm{3}},{\displaystyle }& {\displaystyle }{\mathit{\alpha }}_i(x,z)=\mathrm{0}\text{ otherwise},\end{array}$$

with h₃=1500 m. This initial mixing ratio distribution has been designed to be 𝒞² (and in fact 𝒞³) 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 second-order 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 chemistry-transport studies in academic cases. So far, ToyCTM relies on classical numpy arrays. Its object-oriented 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 Results

3.1 Case 1: thin layer formation under wind shear

Figure 1 Final state of the numerical simulation for Case 1 after simulations. Panels (a–d) show the results obtained with Upwind, VL, PPM and DL99, respectively, and (e) represents the exact solution discretized on model grid. In (a–e), the contour of the exact solution is materialized by a white parallelogram.

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 third-order PPM scheme performing better than the second-order Van Leer scheme and the first-order Upwind scheme in terms of all the diagnostics that have been calculated. More surprisingly, the first-order 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.

Table 3 Performance of simulations performed with the Upwind, VL, PPM and DL99 vertical advection schemes relative to the discretized exact solution for Case 1: percent relative error in $\Vert \cdot {\Vert }_{\mathrm{1}}$ and $\Vert \cdot {\Vert }_{\mathrm{2}}$ and percent of total tracer mass contained in the correct envelope.

3.2 Case 2: long-range 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).

Figure 2 Final state of the numerical simulation for Case 2 after simulations. Panels (a–d) show the results obtained with Upwind, VL, PPM and DL99, respectively; (e) represents the exact solution (which is strictly equal to the initial state due to periodicity in time of the case); and (f) represents the state of the DL99 simulation after 36 h simulation.

Table 4 Performance of simulations performed with the Upwind, VL, PPM and DL99 vertical advection schemes relative to the exact solution for Case 2: percent relative error in $\Vert \cdot {\Vert }_{\mathrm{1}}$ and $\Vert \cdot {\Vert }_{\mathrm{2}}$ and percent of total tracer mass contained in the correct envelope.

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 Convergence rate results for the four tested vertical advection schemes as a function of the number of horizontal points n_x in the horizontal direction, in $\Vert \cdot {\Vert }_{\mathrm{1}}$ (a) and $\Vert \cdot {\Vert }_{\mathrm{2}}$ (b) for Case 3. The black dashed lines show the slope expected for order 1 (long dash) and order 2 (short dash).

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 higher-order 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, higher-order 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 order-1 rate (Fig. A1). Examination of the simulation outputs for these configurations reveal that progress in accuracy is hindered by undesirable small-scale oscillations that degrade accuracy (not shown).

Table 5 Convergence rates in $\Vert \cdot {\Vert }_{\mathrm{1}}$ and $\Vert \cdot {\Vert }_{\mathrm{2}}$.

4 Discussion

The numerical experiments that have been presented confirm the interest of using the Després and Lagoutière (1999) scheme for vertical transport in chemistry-transport 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 long-range 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 chemistry-transport models.

Improved accuracy of a low-order scheme compared to higher-order schemes for a given resolution is not impossible from a theoretical point of view but still counterintuitive since higher-order schemes are designed to reduce numerical error at any given resolution compared to lower-order schemes due to “smarter” reconstruction procedures. Theory imposes that, if model resolution is fine enough and if the tracer field is smooth, higher-order schemes should be more accurate than lower-order schemes. However, as shown by Godunov and Bohachevsky (1959), linear higher-order 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 non-linear slope-limiters 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 slope-limiters 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 higher-order 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 (𝒞³) 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 order-2 schemes if model vertical resolution Δz is equal to 1000 or 500 m, but that due to their faster convergence rate, order-2 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 first-order 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 chemistry-transport models for the advection of polluted plumes thinner than ≃6 Δz (Δz being the model's vertical resolution), while higher-order 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 order-2 schemes in realistic wind conditions including wind shear.

5 Conclusions

The first-order, antidiffusive advection scheme of Després and Lagoutière (1999) has been compared to classical second-order schemes of Van Leer (1977) and Colella and Woodward (1984) in three idealized test cases representing long-range 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 long-range 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, higher-order schemes perform better than the Després and Lagoutière (1999) first-order 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, higher-order 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 chemistry-transport 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 chemistry-transport 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 low-diffusion 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 non-linear 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?

Appendix A: Case 4: convergence test in a smooth case

A1 Case definition

Figure A1 Convergence rate results for the four tested vertical advection schemes as a function of the number of horizontal points n_x in the horizontal direction, in $\Vert \cdot {\Vert }_{\mathrm{1}}$ (a) and $\Vert \cdot {\Vert }_{\mathrm{2}}$ (b) for Case 4.

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:

$$\begin{array}{ll}{\displaystyle }{\mathit{\alpha }}_i(x,z)=& {\displaystyle }{\displaystyle \frac{{\mathit{\alpha }}_m}{\mathrm{4}}}\left(\mathrm{1}+\text{cos}\left({\displaystyle \frac{\mathit{\pi }(z-H/\mathrm{2})}{h_{\mathrm{4}}}}\right)\right)\\ \text{(A1)}& {\displaystyle }& {\displaystyle }\times \left(\mathrm{1}+\text{cos}\left({\displaystyle \frac{\mathit{\pi }(x-L/\mathrm{2})}{\mathit{\delta }}}\right)\right),\end{array}$$

with $h_{\mathrm{4}}=\frac{H}{\mathrm{2}}=\mathrm{6000}\mathrm{m}$ and $\mathit{\delta }=\frac{L}{\mathrm{2}}=\mathrm{5}\times {\mathrm{10}}^{\mathrm{5}}\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):

Table A1 Convergence rates in $\Vert \cdot {\Vert }_{\mathrm{1}}$ and $\Vert \cdot {\Vert }_{\mathrm{2}}$ for Case 4. Convergence rates for the DL99 scheme, marked with a ^* symbol, are evaluated between nx=40 and nx=80, before numerical oscillations appear and begin to degrade the result.