the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
An improved version of the piecewise parabolic method advection scheme: description and performance assessment in a bidimensional test case with stiff chemistry in toyCTM v1.0.1
Sylvain Mailler
Romain Pennel
Laurent Menut
Arineh Cholakian
This study presents a novel method to estimate the performance of advection schemes in numerical experiments along with a semirealistic, nonlinear, stiff chemical system. This method is based on the examination of the “signature function”, an invariant of the advection equation. Apart from exposing this concept in a particular numerical test case, we show that a new numerical scheme based on a combination of the piecewise parabolic method (PPM) with the flux adjustments of Walcek outperforms both the PPM and the Walcek schemes for inert tracer advection as well as for advection of chemically active species. From a fundamental point of view, we think that our evaluation method, based on the invariance of the signature function under the effect of advection, offers a new way to evaluate objectively the performance of advection schemes in the presence of active chemistry. More immediately, we show that the new PPM + W (“piecewise parabolic method + Walcek”) advection scheme offers chemistrytransport modellers an alternative, highperformance scheme designed for Cartesiangrid Eulerian chemistrytransport models, with improved performance over the classical PPM scheme. The computational cost of PPM + W is not higher than that of PPM. With improved accuracy and controlled computational cost, this new scheme may find applications in other fields such as ocean models or atmospheric circulation models.
 Article
(4071 KB)  Fulltext XML
 BibTeX
 EndNote
Chemistrytransport models are models that aim at representing the concentration of trace gases and particles in the atmosphere. Many such tools exist and are used for several purposes, including research and operational forecast. The core of such models consists of a chemical solver adapted to stiff ordinary differential equation (ODE) systems along with a framework for solving the advection equation for all the chemical species.
Among the possible strategies to solve the advection equation in chemistrytransport models are the fluxbased advection schemes, based on the ideas of Godunov (1959), including the Van Leer (1977) and Colella and Woodward (1984) schemes. The Walcek (2000) scheme, an improvement of Van Leer (1977), has also been of common use in chemistrytransport models. For example, GEOSChem provides an advection framework based on the FV3 module implementing the Putman and Lin (2007) method (Martin et al., 2022), based on the Colella and Woodward (1984) scheme for 1D advection. The Colella and Woodward (1984) PPM scheme is also implemented in the CMAQ model (Byun and Schere, 2006; Zhao et al., 2020). The CHIMERE model (Menut et al., 2021) also provides the Van Leer (1977) and Colella and Woodward (1984) schemes for horizontal advection, while vertical advection can be treated with either the Van Leer (1977) scheme or the Després and Lagoutière (1999) antidiffusive advection scheme (Lachatre et al., 2020). The Walcek (2000) advection scheme, an improved version of Van Leer (1977) with reduced numerical diffusion, is also used in chemistrytransport models including CCATTBRAMS (Freitas et al., 2012) and LOTOSEUROS (Timmermans et al., 2022). Therefore, while not the only ones used in chemistrytransport models, the schemes we study here are the base of the numerical resolution of advection in some of the most common chemistrytransport models. Other popular schemes include versions of the Bott (1989) scheme, which is less diffusive than PPM but has the inconvenience of being nonmonotonic and therefore tends to generate extreme values or oscillation in the presence of large concentration gradients (Byun and Schere, 2006), while, by construction, PPM, Van Leer and Walcek schemes enforce mass conservation and monotonicity (Van Leer, 1977; Colella and Woodward, 1984; Walcek, 2000). Despite the availability of higherorder and less diffusive schemes such as the Prather scheme (Prather, 1986) or MPDATA (e.g. Waruszewski et al., 2018), PPM is still considered “highly accurate and efficient enough to be useful” (Harris et al., 2021). A new implementation of the PPM scheme in the CAMx chemistrytransport model has been designed recently to take advantage of GPU computational facilities (Cao et al., 2023). Since the PPM scheme and, to a lesser extent, the Van Leer (1977) and Walcek (2000) schemes are widely used at least in chemistrytransport modelling, it is important to look for ways to improve these schemes while maintaining their desirable properties of robustness and numerical efficiency. In this direction, the goal of the present study is to compare and assess the performance of these schemes in a bidimensional, academic framework including active chemistry; to build an improved version of the PPM scheme, the PPM + W scheme; and to compare the performance of this new scheme to the abovecited classical schemes.
In the past, many studies have focused on developing, improving and evaluating advection schemes (e.g. LeVeque, 1996; Nair and Lauritzen, 2010; Lauritzen et al., 2012, 2014), but very few studies have tackled the evaluation of numerical systems combining advection and chemistry in an academic framework. The most significant step in this direction is the study of Lauritzen et al. (2015), who introduced a toy chemistry scheme mimicking the photolysis and recombination of a virtual stratospheric species. With this simplified nonlinear chemical system, they have tested chemistryadvection combinations, and the errors they generate as diagnostics of massconservation issues. They have noted that combining such a chemical system with the advection solver may reveal problems that are not generated by inert tracer advection. More recently, Lachatre et al. (2022) have shown that changes in the advection formulation may have significant effects on the behaviour of nonlinear chemical processes in the troposphere (in their case, the oxidation pathways of SO_{2} in a midtropospheric volcanic plume). Therefore, like Lauritzen et al. (2015), we feel that it is important to test advection schemes not only with inert tracers but also with active chemistry.
Our goal is to provide such a test case for conditions more representative of tropospheric chemistry at the scale of an urban area and to deploy new tools to evaluate advection schemes in the presence of active chemistry. To meet this objective, apart from classical methods and metrics, we introduce a novel idea, the “signature function”, that permits giving a lower bound of model error compared to the exact solution for problems with inert tracer advection and isolating the error due to advection itself in problems including active chemistry. Even though this method is related to the area coordinate introduced by Nakamura (1996), we introduce a new formulation of this idea along with a way to use it to construct a new error estimate which can be used in problems of pure advection as well as in advection with active chemistry. Apart from this novel way of evaluating advection error in cases with active chemistry, we also propose a new “hybrid” advection scheme, the piecewise parabolic method + Walcek (PPM + W) scheme, made of the PPM scheme with the Walcek (2000) flux adjustments in the vicinity of the extrema.
The PPM + W advection scheme and the concept of the signature function are tested within the toyCTM academic chemistrytransport model, already used in Mailler et al. (2021) to test the use of the antidiffusive scheme of Després and Lagoutière (1999) for vertical advection in the atmosphere. The model version used for this study, toyCTM v1.0 (Mailler and Pennel, 2023), includes horizontal advection with the following schemes (according to user's choice): Godunov (1959), Van Leer (1977), Walcek (2000), Colella and Woodward (1984), and PPM + W (present study); the Després and Lagoutière (1999) scheme is also available for the vertical direction. Chemical processes are solved using an Euler backward iterative (EBI) method (Hertel et al., 1993). As reviewed in Cariolle et al. (2017), this EBI scheme or closely related schemes are used in the MOZART model (Emmons et al., 2010), the ECHAM5HAMMOZ model (Pozzoli et al., 2008), the TM5 model (Huijnen et al., 2010) and the UKCA climatecomposition model (O'Connor et al., 2014; Esentürk et al., 2018).
In Sect. 2, we describe the flux and chemistry of the numerical experiment we have implemented. In Sect. 3, we present the set of simulations we have performed and analysed as well as a description of the advection schemes (Sect. 3.1), the chemical solver (Sect. 3.2) and the timestepping strategy (Sect. 3.3) implemented in toyCTM v1.0. In Sect. 4, we present the concept of signature function that we introduce in this study for the analysis of simulation results. Section 5 compares and discusses the results obtained with the various numerical schemes, and our conclusions are presented in Sect. 6.
2.1 Chemical mechanism
The chemical mechanism used here (Reactions R1–R12) includes a subset of the main reactions of tropospheric gasphase chemistry. Reactions (R1)–(R3) are the three reactions that constitute the Leighton system, and Reactions (R4) to (R7) account for the formation of the hydroxyl radical OH through the photolysis of ozone in the presence of water vapour. Reactions (R8) and (R9) account for the production of the hydroperoxyl radical through oxidation of CO and the oxidation of NO into NO_{2} by HO_{2}. Reactions (R10)–(R11) are “termination reactions” that consume the radical species, and Reaction (R12) describes the final consumption of the NO_{x} species by formation of nitric acid.
The reaction constants of Reactions (R1)–(R12) have been taken mostly from Seinfeld and Pandis (1997), with a temperature of 298 K and pressure of 101 325 Pa. The photolysis rates have been set to typical midday values (e.g. Mailler et al., 2016):

${j}_{\mathrm{R}\mathrm{1}}=\mathrm{8}\times {\mathrm{10}}^{\mathrm{3}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{\mathrm{1}}$

${j}_{\mathrm{R}\mathrm{4}}=\mathrm{2.5}\times {\mathrm{10}}^{\mathrm{5}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{\mathrm{1}}$.
Apart from the chemically active species defined in Reactions (R1)–(R12), we introduce an inert tracer denoted TRC, which undergoes no chemical reaction and is passively advected by the flow. In chemistrytransport models, Reactions (R4)–(R7) are typically lumped in one single reaction, ${\mathrm{O}}_{\mathrm{3}}+h\mathit{\nu}\to \mathrm{OH}+\mathrm{OH}$, with a pseudoreaction rate constant that depends on the concentration of air molecules and of water vapour molecules and on applying the quasisteadystate approximation to O(^{1}D). For this study, we have chosen to treat O(^{1}D) as a prognostic species to preserve the full chemical stiffness of the problem, with lifetimes ranging from $\simeq \mathrm{4}\times {\mathrm{10}}^{\mathrm{9}}$ s for O(^{1}D) to several days for CO.
Of course key processes like oxidation of methane and of other volatile organic compounds are not taken into account in the above mechanism, but it retains some key features of tropospheric chemistry, which we think important:

extreme stiffness;

OH production, which depends on the presence of ozone, water vapour and sunlight;

nonlinear behaviour of ozone production (in this simplified system, ozone production depends on the simultaneous presence of nitrogen oxides, OH and available CO for oxidation).
2.2 Definition of test case
2.2.1 Simulation domain
The simulations are performed on a domain $\mathcal{D}=\left[\mathrm{0},L\right]\times \left[\mathrm{0},L\right]\times \left[\mathrm{0},H\right]$, where L=10^{5} m and H=1000 m. Since we will use only barotropic winds, the problem is in fact bidimensional in x–y, with no z dependance. However, the choice has been made to formally treat the problem as tridimensional in order to be able to use quantities such as air density and reaction rate constants with their usual magnitudes and units. Due to the barotropic nature of the problem, discretization in the vertical direction is in one single cell, while the x and y dimensions are split evenly into n=25 subintervals each. This corresponds to a resolution $\mathit{\delta}x=\mathrm{4}\times {\mathrm{10}}^{\mathrm{3}}$ m, rather typical of regionalscale chemistrytransport modelling. Domain 𝒟 is therefore discretized into n^{2} cells, each cell with thickness H and horizontal section $\mathit{\delta}{x}^{\mathrm{2}}=\frac{{L}^{\mathrm{2}}}{{n}^{\mathrm{2}}}$.
2.2.2 Wind field
The flow we use in this study is the swirling deformational flow introduced by LeVeque (1996) (their Eqs. 9.5–9.6):
with
T is the half period of the experiment, and the design of the flow is such that all fluid particles are back at their original location after time T, but in between they have undergone a deformation, which is maximal at time $\frac{T}{\mathrm{2}}$. Here, while the LeVeque (1996) study formulates the problem with nondimensional scales for time and space, we set a dimensional scale length L=10^{5} m and half period T=86 400 s. The velocity field corresponding to these values is depicted in Fig. 1. Equation (1) ensures that the wind is zero at domain boundaries ($x\in \left\{\mathrm{0},L\right\}$ or $y\in \left\{\mathrm{0},L\right\}$) so that no mass enters or leaves domain 𝒟. Therefore, no boundary conditions for concentrations are needed.
The timedependent streamfunction for this flow is
2.3 Initial conditions
The numerical experiments will be conducted in the domain 𝒟 defined above, with the chemical scheme described above. To define the initial conditions, we introduce a concentration profile (between 0 and 1) as follows:
The initial conditions are defined as follows (in terms of the mixing ratio):

${\mathit{\alpha}}_{\mathrm{TRC}}=\mathrm{100}\phantom{\rule{0.125em}{0ex}}\mathrm{ppb}\times \mathit{\phi}\left(x,y\right)$ (see Fig. 2)

${\mathit{\alpha}}_{{\mathrm{O}}_{\mathrm{3}}}=\mathrm{30}$ ppb

α_{CO}=500 ppb

${\mathit{\alpha}}_{\mathrm{NO}}=\mathrm{100}\phantom{\rule{0.125em}{0ex}}\mathrm{ppb}\times \mathit{\phi}\left(x,y\right)$

${\mathit{\alpha}}_{{\mathrm{NO}}_{\mathrm{2}}}=\mathrm{10}\phantom{\rule{0.125em}{0ex}}\mathrm{ppb}\times \mathit{\phi}\left(x,y\right)$

${\mathit{\alpha}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}=\mathrm{8.044}\times {\mathrm{10}}^{\mathrm{3}}$.
For water vapour, the mixing ratio ${\mathit{\alpha}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}=\mathrm{8.044}\times {\mathrm{10}}^{\mathrm{3}}$ corresponds to a specific humidity of 5 g kg^{−1}. The initial concentrations of the other active species (O, O(^{1}D), OH, HO_{2}, CO_{2}, H_{2}O_{2}, HNO_{3}) are initialized to zero. They will be produced by Reactions (R1), (R4), (R5), (R8) and (R12).
Finally, another species of inert tracer, TRCb, is introduced so that
This species is designed so that, at the initial time and later along the run, in the exact solution (but not necessarily in numerical solutions), all along the run, the sum ${\mathit{\alpha}}_{\mathrm{TRCb}}+{\mathit{\alpha}}_{\mathrm{NO}}+{\mathit{\alpha}}_{{\mathrm{NO}}_{\mathrm{2}}}+{\mathit{\alpha}}_{{\mathrm{HNO}}_{\mathrm{3}}}$ is uniform, constant and equal to 110 ppb.
3.1 Advection schemes
3.1.1 Existing advection schemes
The following existing advection schemes have been tested in the study:
These schemes are fluxbased, upwindbiased, semiLagrangian schemes based on polynomial reconstructions of the average concentrations. These polynomial reconstructions are piecewiseconstant for Godunov (1959), piecewiselinear for Van Leer (1977) and Walcek (2000) (Fig. 3a, b), and piecewise parabolic except in the vicinity of the extrema for PPM (Fig. 3c). The Van Leer (1977) scheme exhibits a discontinuity in the vicinity of the maximum, with the concentrations having a positive jump towards the maximum (Fig. 3a). As a consequence of this discontinuity, due to the upwindbiased strategy, the fluxes going out of the maximum (from the high side of the discontinuity) will be systematically overestimated compared to the fluxes going into the maximum (from the low side of the discontinuity), thereby tending to advect too much mass out of the maximum and not enough mass into the maximum. Walcek (2000) presents his scheme as a way to counteract this bias, by adjusting the flux estimates in the cells next to the maximum, in order to intentionally overestimate the fluxes going into the maximum to counteract the excessive estimation of the fluxes out of the maxima (Fig. 3b).
3.1.2 The PPM + W scheme
The PPM scheme presents the same caveat as Van Leer in the vicinity of extrema, with a strong discontinuity on each side of the extremum (Fig. 4a), with the effect of underestimating the mass flux into the maximum and overestimating the mass flux out of the maximum. Therefore, since Walcek (2000) has proven that his flux adjustments in the vicinity of the extrema are successful in improving the Van Leer scheme, it makes sense to try applying the same flux adjustments to PPM, which seems to have a behaviour similar to Van Leer (1977) in the vicinity of mixingratio extrema (Fig. 3a, c). To test this idea, we design a new scheme based on PPM but applying the Walcek flux adjustments in the vicinity of the extrema (Figs. 3d and 4b). We call this scheme PPM + W, standing for the “piecewise parabolic method + Walcek flux adjustments”. The PPM + W has the same behaviour as Walcek in the extrema and the neighbouring cells and the same behaviour as PPM in all other cells.
We detail here the procedure applied for this scheme. Let (αi) be the values of the mixing ratio in the model cells numbered by the 1D index i, with a Courant number ν. For the sake of simplicity we assume that δx=1. The objective of this procedure is to calculate the average mixing ratio between x_{i+½} and ${x}_{i+\frac{\mathrm{1}}{\mathrm{2}}}\mathit{\nu}$, ${\stackrel{\mathrm{\u0303}}{\mathit{\alpha}}}_{i+\frac{\mathrm{1}}{\mathrm{2}}}$ (the rest of the implementation of advection from this estimate is detailed in Lachatre et al., 2020).
The procedure is as follows.

If $\left({\mathit{\alpha}}_{i+\mathrm{1}}{\mathit{\alpha}}_{i}\right)\left({\mathit{\alpha}}_{i+\mathrm{2}}{\mathit{\alpha}}_{i+\mathrm{1}}\right)>\mathrm{0}$ and $\left({\mathit{\alpha}}_{i\mathrm{1}}{\mathit{\alpha}}_{i\mathrm{2}}\right)\left({\mathit{\alpha}}_{i}{\mathit{\alpha}}_{i\mathrm{1}}\right)>\mathrm{0}$,
the current cell is not a neighbour of a maximum. We estimate ${\stackrel{\mathrm{\u0303}}{\mathit{\alpha}}}_{i+\frac{\mathrm{1}}{\mathrm{2}}}$ following the piecewise parabolic method procedure described in Colella and Woodward (1984). 
Otherwise we estimate the Walcekadjusted flux as follows:
$s=\mathrm{sign}\left({\mathit{\alpha}}_{i+\mathrm{1}}{\mathit{\alpha}}_{i}\right)\mathrm{min}\left(\frac{\mathrm{1}}{\mathrm{2}}\left{\mathit{\alpha}}_{i+\mathrm{1}}{\mathit{\alpha}}_{i\mathrm{1}}\right,\mathrm{2}\left{\mathit{\alpha}}_{i+\mathrm{1}}{\mathit{\alpha}}_{i}\right,\mathrm{2}\left{\mathit{\alpha}}_{i}{\mathit{\alpha}}_{i\mathrm{1}}\right\right)$, the Van Leer slope;
if $\left({\mathit{\alpha}}_{i+\mathrm{1}}{\mathit{\alpha}}_{i}\right)\left({\mathit{\alpha}}_{i+\mathrm{2}}{\mathit{\alpha}}_{i+\mathrm{1}}\right)\le \mathrm{0}$, $\mathit{\beta}=\mathrm{1.75}\mathrm{0.45}\mathit{\nu}$,
else $\mathit{\beta}=\mathrm{max}\left(\mathrm{1.5},\mathrm{1.2}+\mathrm{0.6}\mathit{\nu}\right)$,
${\stackrel{\mathrm{\u0303}}{\mathit{\alpha}}}_{i+\frac{\mathrm{1}}{\mathrm{2}}}={\mathit{\alpha}}_{i}+\frac{\mathrm{1}}{\mathrm{2}}\left(\mathrm{1}\mathit{\nu}\right)\times \mathit{\beta}\times s$.
The β coefficient (β>1 by construction) is introduced by Walcek (2000) to steepen the Van Leer slopes in the vicinity of the maxima to obtain the desired overestimation of tracer fluxes into the maxima. These steepened slopes are visible in Fig. 3b for the Walcek (2000) scheme and in Figs. 3d and 4b for the PPM + W scheme.
3.1.3 Computational cost of the advection schemes
To evaluate the computational cost of these advection schemes, advection of a 1D vector composed of 2×10^{5} cells has been performed over 520 time steps, corresponding to 1.04×10^{8} calls to the reconstruction routine, plus the update of the mixingratio values at each time step. The calculation time for all these advection schemes is presented in Table 1, showing that the schemes using linear reconstruction (Van Leer and Walcek) are less costly than the schemes using parabolic reconstruction (PPM and PPM + W) due to the simpler calculation. Interestingly, the computation cost of PPM + W is slightly smaller than the cost of PPM, possibly because in the cells neighbouring a concentration extremum the reconstruction is linear in PPM + W instead of parabolic in PPM (Fig. 4).
3.1.4 Convergence properties
A first comparison between the PPM + W advection scheme and the other four tested schemes has been performed in terms of numerical convergence. For this purpose, an inert tracer with a squared cosine bell distribution (initially) has been advected over a 1D periodic domain of unit length, with a constant and uniform speed, with a unit duration. This convergence test has been performed dividing the domain into 10, 20, 40, 80, 160 and 320 grid cells and with a Courant number of 0.5. The results of this convergence test are shown in Fig. 5. Several features can be observed in Fig. 5. The convergence rates, defined as the opposite of the log–log slope between the two last data points (n_{x}=160 and n_{x}=320), are given in Table 2. From these results, several observations can be made.
First, the PPM + W advection scheme performs better on this simple convergence test than the other tested schemes. In particular, throughout all the resolution range, the error obtained in both ∥⋅∥_{1} (Fig. 5a) and ∥⋅∥_{2} (Fig. 5b) with the PPM + W scheme is 30 % to 50 % lower than with the classical PPM scheme. This difference persists even for high resolutions. On the contrary, the Walcek scheme strongly outperforms the Van Leer scheme for coarse resolutions, but this difference tends to diminish for higher resolution. In other terms, it looks like the Walcek flux corrections when applied to the Van Leer scheme permit the improvement of accuracy when model resolution is coarse but that the same flux corrections when applied to the PPM scheme improve accuracy even for fine resolutions. Finally, the convergence rate for PPM + W is similar to that of PPM, in our case around 2.5 for the ∥⋅∥_{1} convergence rate and around 2 for ∥⋅∥_{2} (Fig. 2).
3.2 Chemistry solver: the Euler backward iterative method
The stiff chemical system is integrated using an EBI method. As described in Hertel et al. (1993) and Cariolle et al. (2017), we obtain the concentration vector c(t+δt_{chem}) as the solution of
For the present study, the focus is to test the performance of the advection scheme in articulation with active chemistry. Due to this focus, we limit errors in the resolution of the chemical system using a short time step for chemistry (δt_{chem}=20 s).
Equation (7) is a nonlinear, fixedpoint equation and can be solved only numerically, usually with an iterative method. Formally, Eq. (7) guarantees exact mass conservation. However, this is true only if a very good convergence of the solution is reached (Cariolle et al., 2017). To limit violation of mass conservation in our study, we have set a very strict convergence criterion for the iterative resolution of Eq. (7), stopping iteration when the estimate of c(t+δt_{chem}) yields a relative difference less than ε for each species between c(t+δt_{chem}) and $\frac{\mathit{c}\left(t\right)+\mathit{\delta}{t}_{\mathrm{chem}}\mathit{P}\left(t+\mathit{\delta}{t}_{\mathrm{chem}}\right)}{\mathrm{1}+\mathit{\delta}{t}_{\mathrm{chem}}\mathit{L}\left(t+\mathit{\delta}{t}_{\mathrm{chem}}\right)}$. The convergence parameter is set to $\mathit{\epsilon}={\mathrm{10}}^{\mathrm{4}}$ in the UKCA chemistrytransport model (Esentürk et al., 2018) and $\mathit{\epsilon}={\mathrm{10}}^{\mathrm{6}}$ in the present study. This very strict convergence criterion is in line with the short chemical time step to obtain the bestpossible numerical solution of the chemical evolution of the system, even at the cost of slow computations. Again, this choice is due to the purpose of this study to test the performance of the advection scheme, limiting as much as possible the errors in the chemical solver.
3.3 Time stepping
The advection time step δt_{adv} has been set to δt_{adv}=1800 s. With maximal wind module U≃1.8 m s^{−1} (Fig. 1) and $\mathit{\delta}x=\mathrm{4}\times {\mathrm{10}}^{\mathrm{3}}$ m, this yields a maximal Courant number ν_{max}≃0.8. The timestepping strategy follows Strangstyle time stepping (Strang, 1968), with the steps as follows:

Integrate chemistry over $\frac{\mathit{\delta}{t}_{\mathrm{adv}}}{\mathrm{2}}$ (45 chemical time steps).

Integrate zonal advection over $\frac{\mathit{\delta}{t}_{\mathrm{adv}}}{\mathrm{2}}$.

Integrate meridional advection over δt_{adv}.

Integrate zonal advection over $\frac{\mathit{\delta}{t}_{\mathrm{adv}}}{\mathrm{2}}$.

Integrate chemistry over $\frac{\mathit{\delta}{t}_{\mathrm{adv}}}{\mathrm{2}}$ (45 chemical time steps).
Table 3 summarizes the six simulations that have been performed. The abovedescribed case with Reactions (R1)–(R12) has initial conditions as described in Sect. 2.3, a domain shape and discretization as described in Sect. 2.2.1, advection schemes as described in Sect. 3.1, and a chemical solver as described in Sect. 3.2. Along with the simulations performed with each of the five advection schemes, a “base” simulation has been performed with the same setup as the other simulations but without advection. As discussed later (Sect. 4), this base simulation will serve as a benchmark to estimate advection errors in the other five simulations.
Godunov (1959)Van Leer (1977)Walcek (2000)Colella and Woodward (1984)3.4 Conservation properties
It is worth noting that, by construction, fluxbased advection integration is massconservative since the mass flux out of a cell through a facet is compensated exactly by the mass flux into the neighbouring cell through the same facet. Equation (7) also guarantees mass conservation in the Euler backward scheme as soon as the chemical reactions themselves are balanced (Hertel et al., 1993) (which is the case of Reactions R1–R12 except for the imbalance of dioxygen in Reaction R8 due to integrating reaction $\mathrm{H}+{\mathrm{O}}_{\mathrm{2}}\to {\mathrm{HO}}_{\mathrm{2}}$ into the kinetically limiting time step $\mathrm{CO}+\mathrm{OH}\to {\mathrm{CO}}_{\mathrm{2}}+\mathrm{H}$).
However, due to the finite number of iterations in the iterative resolution of Eq. (7), mass conservation is only enforced with a finite precision of $\mathit{\epsilon}={\mathrm{10}}^{\mathrm{6}}$ (see Sect. 3.2). Therefore, the relative mass imbalance in the outputs for C, active N (without taking into account N_{2}) and H can be expected to be of the order of ε, since the Euler backward iterative scheme (Eq. 7) is in principle massconservative, but mass conservation is obtained “only if a good convergence of the solution is reached” (Cariolle et al., 2017).
Mass calculations for C, active N, H and TRC have been performed between the beginning and the end of the simulations. The results of this calculation for the PPM + W and base simulations are given in Table 4, showing that the relative mass imbalance at the end of simulation is $\simeq {\mathrm{10}}^{\mathrm{6}}$ for active N, $\simeq {\mathrm{10}}^{\mathrm{12}}$ for C and H, and $\simeq {\mathrm{10}}^{\mathrm{15}}$ for TRC (which has no chemistry and is therefore affected only by advection which, as discussed above, ensures mass conservation up to numerical accuracy for the fluxbased schemes we have implemented).
The imbalance results are similar for all simulations, including the base simulation that has no advection, which shows that the small mass imbalance for chemically active species (up to $\simeq {\mathrm{10}}^{\mathrm{6}}$) is essentially due to the finite precision in the Euler backward iterative chemistry solver.
In summary, the integration strategy we introduce above permits the conservation of mass (up to numerical accuracy for inert species and up to an arbitrary numerical tolerance defined by the user, in our case $\mathit{\epsilon}={\mathrm{10}}^{\mathrm{6}}$, for chemically active species) and conserves initially uniform mixing ratios up to numerical accuracy for the species with no active chemistry. All the advection schemes presented above are also built to respect monotonicity: they do not create new mixingratio extrema.
4.1 Accuracy of inert tracer concentrations
As the LeVeque (1996) flow is designed so that at t=T every Lagrangian particle is back at its original location, it is possible to estimate the accuracy of numerical simulation by comparing the final simulated tracer concentration field to its initial value, therefore giving access to the magnitude of numerical error.
In the present study, we will estimate model error in ∥⋅∥_{1}, introducing E_{1} as the normalized $\Vert \cdot {\Vert}_{\mathrm{1}}$ error in the mixing ratio:
where index i spans the entire domain. In the present case where density ρ_{i} and cell volume do not vary across cells, Eq. (8) boils down to
4.2 The signature function for inert tracer advection
Here we introduce a new idea to evaluate advection schemes. As far as we know, this idea has not been tested in the past literature but resembles the areacoordinate formulation used by Nakamura (1996).
Let us imagine a fluid with density $\mathit{\rho}\left(x;y;z;t\right)$ in a threedimensional domain 𝒟, advecting a tracer having initially a mixing ratio $\mathit{\alpha}\left(x;y;z;t=\mathrm{0}\right)\in \left[\mathrm{0},\mathrm{1}\right]$, following equations
For any given time t and any mixing ratio $\mathrm{0}\le X\le \mathrm{1}$, we can define 𝒮^{t}(X) as the mass of fluid in the volume 𝒟^{t}(X) defined as the set of all $\left(x;y;z\right)$ where the tracer mixing ratio $\mathit{\alpha}\left(x;y;z;t\right)<X$, divided by the entire mass of fluid in 𝒟:
where H is the Heaviside step function (H(u)=1 if u>0; H(u)=0 if u≤0). The 𝒮^{t} function can be, in some sense, interpreted as a massweighted cumulative probability density function of the tracer mixing ratio. If we reduce this definition to 2D flows with uniform density, 𝒮^{t} is related to the reciprocal function of the areacoordinate formulation of Nakamura (1996).
With this definition, we always have S^{t}(0)=0 (for all t values) and S^{t}(X)→1 when $X\to +\mathrm{\infty}$ (more precisely, S^{t}(X)=1 as soon as X is larger than the maximum value of $\mathit{\alpha}\left(x;y;z;t=\mathrm{0}\right)$ over domain 𝒟).
Equation (12) makes clear that function 𝒮^{t} is invariant during the motion: for any given value of X, 𝒮^{0}(X) is the (normalized) mass of the fluid parcel 𝒟^{0}(X) that has a tracer mixing ratio α<X at t=0. We can observe that, since mixing ratio in Lagrangian parcels is preserved by pure advection (Eq. 11), 𝒟^{t}(X) and 𝒟^{0}(X) represent the same Lagrangian fluid parcel at a different time. We also know that the mass of fluid in Lagrangian parcels is constant in time due to mass conservation for the carrier fluid (Eq. 10). Since the total mass of fluid in 𝒟 is also constant in time for the same reason, this implies that, for all t and X values, 𝒮^{t}(X)=𝒮^{0}(X). In other words, the signature function 𝒮^{t} is an invariant in time of the advection equation. The invariance of this function holds for both divergent and nondivergent flows.
Therefore, since we know that, for the exact solution, 𝒮^{t}=𝒮^{0}, the departure of the numerical evaluation of function 𝒮^{t} from 𝒮^{0} in numerical simulations can be used as an objective measurement of discretization errors.
In practice, in an Eulerian model discretized into N cells, each cell has an evaluation of the tracer mixing ratio ${\mathit{\alpha}}_{i}^{\mathrm{t}}$. The evaluation of 𝒮^{t} is straightforward:
Generally, the numerical evaluation of 𝒮^{t} in an Eulerian model will differ from 𝒮^{0}: the signature of the initial tracer distribution will evolve under the effect of the errors of the advection scheme, and the magnitude of the signature modification can serve as a measure of the extent of advection error. In the particular case in which the flow is nondivergent and the carrier fluid mass ρ_{i}V_{i} is the same in all model cells, the norm1 difference between 𝒮^{t} and 𝒮^{0} can be calculated as
where ${\stackrel{\mathrm{\u0303}}{{\mathit{\alpha}}^{\mathrm{t}}}}_{i}$ is the vector of all mixing ratios in model cells sorted in increasing order. This simplified formula holds only in the specific case in which the flow is invariant and the mass of fluid is distributed evenly between all the model cells. In this case, and only in this case, it is also convenient to introduce the normalized norm1 difference ℰ_{1} between 𝒮^{t} and 𝒮^{0} as
Figure 6 shows an example of how comparing the 𝒮^{t} signature function with 𝒮^{0} permits comparison of the accuracy of two simulations at a time when no analytic solution is easily accessible. Figure 6a–b show the mixing ratio for TRC at $\frac{T}{\mathrm{2}}$ in the Godunov and Van Leer simulations, respectively. Without access to the exact solution, it is hard to compare quantitatively the quality of these simulations at that stage, even though indicators such as the maximal tracer mixing ratio can give partial information. Figure 6c–d show the ${\mathcal{S}}^{\frac{T}{\mathrm{2}}}$ function compared with 𝒮^{0} for the Godunov and Van Leer simulations, respectively. This graphical representation permits giving an intuitive meaning to the normalized signature error ℰ_{1} as the total area between the representative curves of 𝒮^{t} and 𝒮^{0} (shaded in Fig. 6c–d), divided by the total area left of the representative curve of 𝒮^{0} (hatched in Fig. 6c–d). This measure gives an indication of model error based on not only one particular point representative of a part of the tracer distribution (e.g. the maximal value, a very partial indicator), but also an integrated error indicator taking into account the maximal and minimal values as well as the entire tracer distribution in between these values. In the case presented in Fig. 6c–d, the area between 𝒮^{0} and 𝒮^{t} is smaller in the Van Leer simulation (Fig. 6d, with ℰ_{1}=0.218) than in the Godunov simulation (Fig. 6c, with ℰ_{1}=0.571).
4.3 The signature function for advection of active species
The concentrations of active species evolve not only under the effect of advection but also due to chemical reactions. Therefore, the time invariance of the signature function does not hold for these species. However, the signature function can still be used to compare and evaluate simulations if one remarks that, in the case we study here with no variations in air density and air temperature, and with no emissions, the chemistry that takes place in each Lagrangian air parcel is independent of its position. Therefore, for all species, the signature function at time t, 𝒮^{t}, should theoretically be the same as in the base simulation with no advection. If we note ${\stackrel{\mathrm{\u203e}}{\mathcal{S}}}^{\mathrm{t}}$, the signature function at time t in the simulation without advection, for all t values we have ${\mathcal{S}}^{\mathrm{t}}={\stackrel{\mathrm{\u203e}}{\mathcal{S}}}^{\mathrm{t}}$: at any time, the signature of the distribution of every chemical species should be the same in the case with advection as in the case without. In other words, advection should only deform the map of all chemical species and not change the chemistry within each Lagrangian air parcel. With a nonlinear chemical system such as Reactions (R1)–(R12), there is no easy access to the exact solution of the system even without advection: ${\stackrel{\mathrm{\u203e}}{\mathcal{S}}}^{\mathrm{t}}$ is not known exactly for the chemically active species. However, the base simulation resolves the chemical reactions with exactly the same chemical solver as the other simulations but without motion. Therefore, comparing the distribution of a chemical species at time t in a simulation with advection to its distribution in the base simulation without advection will permit having an estimate of how the numerical errors in advection affect the distribution of chemically active species.
As an illustration of this, Fig. 7 shows the ozone mixing ratio at $\frac{T}{\mathrm{2}}$, as simulated in three simulations with advection (Fig. 7a–c) and the base simulation without advection (Fig. 7d). This figure alone does not make it easy to discriminate between the three numerical simulations. Figure 8 shows the signature function for the same three simulations (Van Leer, PPM and PPM + W): this time, both the visualization of the agreement between 𝒮^{t} and ${\stackrel{\mathrm{\u203e}}{\mathcal{S}}}^{\mathrm{t}}$ and the objective calculation of ℰ_{1} clearly show that the PPM + W simulation agrees better with the base simulation in terms of the signature function (ℰ_{1}=0.0622), followed by the PPM simulation (ℰ_{1}=0.0679) and the Van Leer simulation (ℰ_{1}=0.0787). Apart from this quantitative agreement, some qualitative and local conclusions can be drawn from the graphical comparisons between 𝒮^{t} and ${\stackrel{\mathrm{\u203e}}{\mathcal{S}}}^{\mathrm{t}}$, such as the fact that the representation of the ozone minimum in the PPM + W simulation is more accurate than in the PPM simulation, which is visible in the initial part of the signature function (${\mathit{\alpha}}_{{\mathrm{O}}_{\mathrm{3}}}<\mathrm{30}$ ppb) where the shaded area is smaller in PPM + W (Fig. 7c) than in PPM and Van Leer.
5.1 Results for classical performance indicators
Figure 9 shows the evolution of some performance indicators along the course of the experiments. For the same reasons as for the signature function, these metrics should be the same in all simulations in the absence of numerical errors in the representation of advection, so the differences between the time series obtained in the simulations with advection and the base simulation reveal the effects of numerical errors in the representation of advection. It is interesting to note that all the metrics represented in Fig. 9 can be derived directly from the signature function.
Regarding ozone extrema throughout the simulation (Fig. 9a), we observe smaller differences between the different simulations. The PPM + W, PPM and Walcek simulations produce very comparable values of ozone maxima. Surprisingly, the simple Godunov schemes represents the ozone maximum slightly better than all other schemes except Walcek towards the end of the simulation. However, the representation of the ozone minimum by the Godunov simulation is very bad, failing to represent a minimum lower than the background value towards the end of the simulation. The PPM + W and Walcek simulations represent the smallest (and closest to the base simulation) minimum value for O_{3}, followed by PPM and Van Leer.
Regarding the preservation of tracer maxima (Fig. 9b), the PPM + W simulation performs best, with a clear edge over the Walcek and PPM simulations (Fig. 9b). For this metric, the Walcek simulation gives better results than the PPM simulation, with a smaller half plume, closer to the base simulation, Walcek and PPM. The PPM and Walcek simulations perform similarly for this metric.
Figure 9c shows evolution of the nondimensional half volume of the tracer plume during the simulation (defined as the smallest volume containing half of the mass of TRC divided by the total volume of domain 𝒟). As discussed in Lachatre et al. (2020), this is a measure of the diffusivity of the advection schemes. With this metric, we see that the PPM + W simulation performs best, followed by Walcek and PPM.
5.2 Accuracy and signature error
Unlike the partial metrics presented in Fig. 9, the normalized $\Vert \cdot {\Vert}_{\mathrm{1}}$ signature error in the mixing ratios ℰ_{1} is a performance diagnostic for the simulations that takes into account the entire distribution of mixing ratios and not just particular values such as the maximum or minimum. Figure 10 shows ℰ_{1} for TRC, O_{3}, NO and NO_{2} for all the simulations, with the smallest ℰ_{1} indicating the bestperforming simulation.
The first information we get from Fig. 10 is that, for all these compounds, the PPM + W simulation performs best in this regard. For the case of ozone (Fig. 10b), the ℰ_{1} time series discriminates much more between the different simulations, with a clear superiority of the PPM + W simulation over the Walcek and PPM simulations, while the differences between these three simulations in the ozone min–max plot (Fig. 9b) appeared small. Analysing the evolution of ℰ_{1} throughout the simulation shows that, for TRC, NO and NO_{2}, the bestperforming simulations are PPM + W, Walcek and PPM, in this order, but that for ozone, PPM performs better than Walcek, a conclusion that could not be drawn from the min–max plot (Fig. 9b), which indicated a better performance of Walcek regarding both the ozone maximum and the ozone minimum.
While we have shown so far that the analysis of ℰ_{1} permits the drawing of clearer conclusions regarding the performance of the various simulations, in the present case we can confirm these results by comparing ℰ_{1} to a more classical metric, the normalized $\Vert \cdot {\Vert}_{\mathrm{1}}$ error E_{1} (Eq. 9) of the simulations with advection relative to the base simulation without advection. Unlike ℰ_{1}, E_{1} can only be calculated at the final time step (when all the Lagrangian particles are back at their initial location), since the exact solution for t<T is not accessible. Several observations can be made from the values of E_{1} and ℰ_{1} in Table 5.
First, in all cases we always have ℰ_{1}<E_{1}. Qualitatively, this can be interpreted as ℰ_{1} being a weakened form of the $\Vert \cdot {\Vert}_{\mathrm{1}}$ error, retaining the differences in the distribution of mixing ratios but eliminating the differences in the location of the tracer plume.
Interestingly, the performance ranking between the five simulations obtained by analysing the signature error ℰ_{1} is the same for all variables as with E_{1}: for TRC, NO and NO_{2}, PPM + W performs best, followed in this order by Walcek, PPM, Van Leer and Godunov, but for O_{3}, PPM performs better than Walcek. In our study case, analysing ℰ_{1} permits comparison of the performance of the various simulations without access to the exact solution and gives the same results as the analysis of E_{1} (which requires access to the exact solution). Even for simulations with comparable performance to the Walcek and PPM simulations, the signature error ℰ_{1} permits diagnosing which of the two simulations performs better for each variable. Even though the differences between these two simulations are not drastic and depend on the species of interest, the conclusions drawn from the analysis of the signature error ℰ_{1} are the same and correspond to E_{1}. This gives confidence in the usability of the signature error ℰ_{1} as a proxy for simulation accuracy when the exact solution is not available.
Having verified this, it is useful to go back to Fig. 10 and interpret the evolution of ℰ_{1} in time as a hint of when error appears along the simulation. In this regard, we see two very different behaviours between the analysed variables. For TRC, NO and NO_{2}, substantial errors appear almost instantly after 1800 s of simulation (one single time step). This is due to the action of the wind field on the initially very sharp peak of these species (Fig. 2). On the contrary, for O_{3}, having an initially uniform distribution, errors due to advection appear only when sufficient heterogeneity is introduced into the O_{3} map by chemical processes, since all the advection schemes are built to advect exactly a uniform mixing ratio, maintaining its uniformity. Therefore, the onset of advection errors in O_{3} is much slower than in the three species that have initially heterogeneous distributions.
5.3 Results at 1 km resolution
To test the impact of higher resolution on our results, we have performed the same test case as above but refining the resolution from 4 to 1 km and accordingly reducing the time step from 1800 to 450 s. The results for this higherresolution simulation are shown in Table 6. This table shows that, also at this resolution, the PPM + W offers the best performance of all the tested schemes, but, unlike in the 4 km resolution case, this does not hold for all species: at 1 km resolution, the PPM scheme performs comparably to PPM + W for the inert tracer TRC and slightly better for NO_{2}. However, PPM + W clearly performs better for O_{3} and NO. These results suggest that the improvement brought by PPM + W over PPM may tend to become smaller for higher resolutions (contrary to what we observe in Fig. 5).
5.4 Conservation of sum species
As shown by Godunov (1959), “among [linear] schemes of second order accuracy [for the advection equation], there is none which satisfies the monotonicity condition”. As a consequence of this result, known as the Godunov theorem, secondorder, monotonic schemes for the advection equation are all nonlinear. In particular, among all the schemes we study here, all are nonlinear due to their specific treatment of the maxima to preserve monotonicity. However, these nonlinear effects should be kept as small as possible. To test this property, we have tested the conservation of two sum species, which should stay uniform and constant along the run:

The first is ${\mathit{\alpha}}_{\mathrm{TRCb}}+{\mathit{\alpha}}_{\mathrm{NO}}+{\mathit{\alpha}}_{{\mathrm{NO}}_{\mathrm{2}}}+{\mathit{\alpha}}_{{\mathrm{HNO}}_{\mathrm{3}}}$, where TRCb is an inert tracer whose initial distribution is defined by Eq. (6). The mixing ratio of this sum species should stay constant, uniform and equal to 110 ppb.

The second is ${\mathit{\alpha}}_{\mathrm{CO}}+{\mathit{\alpha}}_{{\mathrm{CO}}_{\mathrm{2}}}$, which should stay constant, uniform and equal to 500 ppb due to the conservation of carbon.
Figure 11 shows that, as imposed by the Godunov theorem, only the Godunov (1959) scheme preserves exactly the uniformity of sum species. Other schemes are higherorder and monotonic, but therefore nonlinear, and do not preserve the uniformity of these sum species. The magnitude of the departure of the mixing ratio for these species from its theoretical value is a measure of the nonlinearity of advection schemes, which is an undesirable property since the advection equation itself is linear. Figure 11 shows that nonlinearities of up to 1 %–2 % appear for all schemes in the representation of CO+CO_{2} and that they tend to be strongest in the Walcek and the PPM + W schemes. Regarding $\mathrm{TRCb}+\mathrm{NO}+{\mathrm{NO}}_{\mathrm{2}}+{\mathrm{HNO}}_{\mathrm{3}}$, the nonlinearities reach up to 10 % of the expected value and are strongest in the PPM + W scheme. These results suggest that nonconservation of sum species in the vicinity of mixingratio extrema, which is observed in all secondorder monotonic schemes, tends to be stronger when the Walcek (2000) flux corrections are applied, which might represent a drawback of PPM + W compared to PPM. In spite of (or due to) this increasingly nonlinear behaviour in the vicinity of the extrema, PPM + W preserves the values of extrema better than PPM (Fig. 9b).
We have introduced the signature function 𝒮^{t} as a sort of massweighted cumulative probability distribution function of the tracer mixing ratio (Eq. 12) and shown that 𝒮^{t} is an invariant of the advection equation. This invariant is not a scalar like other classical invariants (tracer mass, minimum and maximum mixing ratio, etc.) but a function and therefore contains much more information than the abovecited, more classical invariants. In fact, these invariants (tracer mass, minimum and maximum mixing ratios, etc.) can be derived directly from 𝒮^{t} so that 𝒮^{t} can be considered a “stronger” invariant. Like other invariants such as the minimum and maximum values of the mixing ratio, it is usually not preserved perfectly by the advection schemes, and the degree of nonconservation of 𝒮^{t} gives a proxy for model error. Based on this idea, we propose the normalized $\Vert \cdot {\Vert}_{\mathrm{1}}$ signature error ℰ_{1} (Eq. 15) as a measure of model error. Graphically, ℰ_{1} is the normalized area between the simulated curve of 𝒮^{t} and its theoretical curve (e.g. Fig. 6c–d).
In this context, we have shown that the signature function and its normalized $\Vert \cdot {\Vert}_{\mathrm{1}}$ error ℰ_{1} can also be used as an error estimate for chemically active species, even though 𝒮^{t} is not an invariant of the system in this case. This is achieved by comparing a simulation including chemistry and advection to a “companion” simulation with chemistry only but no advection (the base simulation in the present study). For each chemical species and at any time t, these two simulations should in theory have the same 𝒮^{t}, but in practice they do not due to errors in advection. Therefore, errors due to advection can be estimated even for chemically active species by calculating the ℰ_{1} error between the two functions (Fig. 8).
We have used this new invariant in order to evaluate a new advection scheme that we have designed for this study, based on the PPM scheme (Colella and Woodward, 1984) with flux corrections in the vicinity of the extrema as in Walcek (2000) (Figs. 3–4). This new advection scheme, which we propose to name PPM + W (for piecewise parabolic method + Walcek), has been tested for both an inert tracer and chemically active species, with a velocity flux from LeVeque (1996) (Eq. 1). A simplified chemical scheme with 12 reactions has been designed, representing daytime photochemistry of nitrogen oxides in the presence of CO, including shortlived species such as O(^{1}D) or O (Reactions R1–R12). With an initial peak of NO and NO_{2} concentrated over a small area (Fig. 2) and initially uniform concentrations of CO and O_{3}, this case generates a sharp ozone minimum colocalized with the NO_{x} peak, a large area with background O_{3} concentration and a belt of high ozone concentrations in between (Fig. 7). With this test case, we have evaluated the PPM + W scheme along with the PPM scheme, the Walcek (2000) scheme, the Van Leer (1977) scheme and the Godunov (1959) scheme. In this case, we have shown that, for all species and all the metrics we have tested, the PPM + W scheme performs better than all the other tested schemes (Table 5), with a normalized $\Vert \cdot {\Vert}_{\mathrm{1}}$ error E_{1} of 16.4 % in O_{3} (19.3 % for PPM), 20.7 % in TRC (24.3 % for Walcek), 28.3 % in NO_{2} (31.2 % for Walcek) and 37.2 % in NO (42.5 % for Walcek). Table 5 also shows the ℰ_{1} signature error, which is always smaller than E_{1}. Interestingly, examining the ℰ_{1} errors for the same variables and the same simulations yields exactly the same conclusions as examining the normalized $\Vert \cdot {\Vert}_{\mathrm{1}}$ error E_{1}. This shows that, even without access to an exact solution, the ℰ_{1} signature error permits comparison of the simulations against each other for inert and reactive species, giving the same conclusions as an accuracy analysis with E_{1}. This being shown, Fig. 10 permits visualization of the evolution of error along the course of the simulation, while E_{1} can be calculated only at t=T because the LeVeque (1996) flux field is designed to guarantee that at that time all the Lagrangian particles should be back at their initial locations.
Therefore, the conclusion of this study is twofold. First, regarding the signature function as an invariant of the advection equation, we feel that this invariant contains much more information than other invariants that have been typically used to check advection schemes, such as the minimum or maximum values of mixing ratios, while not requiring access to the exact solution. In the case of chemistrytransport models, generalizing this concept to more dimensions (by studying the massweighted probability distribution function of all species simultaneously rather than one signature function per species) could be promising. For the same reasons as exposed above, this multidimensional probability distribution function should be an invariant of the advection equation. This concept could be explored to quantify model error in a synthetic way across all variables, instead of separately for each variable. The approach introduced here with the ℰ_{1} signature error could possibly be generalized using statistical tools such as the Wasserstein distance to compare these multidimensional probability distribution functions with each other. The main limit to the use of the signature function as a tool to evaluate the advection framework in realworld geophysical models is that it relies on mixingratio conservation, which holds for pure advection but does not hold in the presence of mixing or diffusion. However, some geophysical compartments like the deep ocean or the stratosphere are substantially affected by mixing only for very long timescales, so the signature function could be a useful tool to verify the behaviour of advection frameworks in such compartments.
From a more applied point of view, the PPM + W advection scheme introduced here performs better than both the Walcek (2000) scheme and the PPM scheme, with all the metrics we have tested and for both inert and active species. It not only preserves the tracer maxima better than the Walcek (2000) scheme (Fig. 9), but also is more accurate than the PPM scheme for the representation of ozone and other active species (Table 5). Of course, these results are proven in the present study only on one 2D test case with active chemistry and on a 1D convergence test. While the results of both the test case and the convergence test consistently indicate a better performance of PPM + W compared to the other advection schemes tested here, they do not include the full range of Courant numbers and tracer patterns that occur in realistic models. For example, an additional numerical experiment presented here (Sect. 5.3 and Table 6) shows that, when the resolution is refined four times compared to our main test case, the advantage of PPM + W over PPM seems to be reduced and, for one species (NO_{2}), PPM + W is even outperformed by PPM in terms of accuracy. However, even in this highresolution test case, the performance of PPM + W is better than that of PPM for the three other tested species. The convergence test performed in this study (Fig. 5) suggests that there is no systematic reduction in the performance of PPM + W at higher resolution. A possible drawback of the PPM + W scheme when compared to the PPM scheme is its stronger nonlinearity in the vicinity of maxima (Fig. 11), which permits better conservation of the maxima themselves but induces more numerical artefacts in the representation of sum species in the vicinity of extrema.
We have also observed (Table 1) that the computation cost of PPM + W is slightly lower than the cost of the PPM scheme, which is used in some of the most popular chemistrytransport models. The improved performance in terms of accuracy and of preservation of tracer extrema without increasing the computational load makes this scheme a very interesting option for chemistrytransport models, in an effort to reduce numerical diffusion, which is important in particular in the presence of nonlinear chemistry as discussed in Lachatre et al. (2022).
This study has been performed using toyCTM v1.0.1 (https://doi.org/10.5281/zenodo.10018706, Mailler and Pennel, 2023). All the Python scripts used to launch the model and to perform the postprocessing of model outputs to obtain the figures in this paper and the numbers in Table 5 are available from https://doi.org/10.5281/zenodo.10018761 (Mailler et al., 2023).
The Fortran code AdvBench v1.0.0 used to evaluate advection performance (Table 1) is available from https://doi.org/10.5281/zenodo.7937121 (Mailler, 2023).
The model toyctm v1.0.1, AdvBench v1.0.0 and all the scripts used to launch the model and postprocess its outputs for the present study are distributed under the GNU General Public License v2.0. No datasets were used for this study.
All the authors contributed to the design of the simulated cases; SM performed and analysed the simulations; SM developed the software with RP.
The contact author has declared that none of the authors has any competing interests.
Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
The numerical calculations in this study benefited from the IPSL Data and Computing Center ESPRI, which is supported by CNRS, SU, CNES and École Polytechnique.
This research has been supported by the Agence de l'Environnement et de la Maîtrise de l'Énergie (ESCALAIR).
This paper was edited by Sylwester Arabas and reviewed by Christopher Walcek, Hilary Weller, and one anonymous referee.
Bott, A.: A Positive Definite Advection Scheme Obtained by Nonlinear Renormalization of the Advective Fluxes, Mon. Weather Rev., 117, 1006–1016, https://doi.org/10.1175/15200493(1989)117<1006:APDASO>2.0.CO;2, 1989. a
Byun, D. and Schere, K. L.: Review of the Governing Equations, Computational Algorithms, and Other Components of the Models3 Community Multiscale Air Quality (CMAQ) Modeling System, Appl. Mech. Rev., 59, 51–77, https://doi.org/10.1115/1.2128636, 2006. a, b
Cao, K., Wu, Q., Wang, L., Wang, N., Cheng, H., Tang, X., Li, D., and Wang, L.: GPUHADVPPM V1.0: a highefficiency parallel GPU design of the piecewise parabolic method (PPM) for horizontal advection in an air quality model (CAMx V6.10), Geosci. Model Dev., 16, 4367–4383, https://doi.org/10.5194/gmd1643672023, 2023. a
Cariolle, D., Moinat, P., Teyssèdre, H., Giraud, L., Josse, B., and Lefèvre, F.: ASIS v1.0: an adaptive solver for the simulation of atmospheric chemistry, Geosci. Model Dev., 10, 1467–1485, https://doi.org/10.5194/gmd1014672017, 2017. a, b, c, d
Colella, P. and Woodward, P. R.: The piecewise parabolic method (PPM) for gasdynamical simulations, J. Comput. Phys., 11, 38–39, https://doi.org/10.1016/00219991(84)901438, 1984. a, b, c, d, e, f, g, h, i, j
Després, B. and Lagoutière, F.: Un schéma non linéaire antidissipatif pour l'équation d'advection linéaire, CR Acad. Sci. I – Math., 328, 939–943, https://doi.org/10.1016/S07644442(99)803012, 1999. a, b, c
Emmons, L. K., Walters, S., Hess, P. G., Lamarque, J.F., Pfister, G. G., Fillmore, D., Granier, C., Guenther, A., Kinnison, D., Laepple, T., Orlando, J., Tie, X., Tyndall, G., Wiedinmyer, C., Baughcum, S. L., and Kloster, S.: Description and evaluation of the Model for Ozone and Related chemical Tracers, version 4 (MOZART4), Geosci. Model Dev., 3, 43–67, https://doi.org/10.5194/gmd3432010, 2010. a
Esentürk, E., Abraham, N. L., ArcherNicholls, S., Mitsakou, C., Griffiths, P., Archibald, A., and Pyle, J.: QuasiNewton methods for atmospheric chemistry simulations: implementation in UKCA UM vn10.8, Geosci. Model Dev., 11, 3089–3108, https://doi.org/10.5194/gmd1130892018, 2018. a, b
Freitas, S. R., Rodrigues, L. F., Longo, K. M., and Panetta, J.: Impact of a monotonic advection scheme with low numerical diffusion on transport modeling of emissions from biomass burning, J. Adv. Model. Earth Sy., 4, M01001, https://doi.org/10.1029/2011MS000084, 2012. a
Godunov, S. K.: Finite difference method for numerical computation of discontinuous solutions of the equations of fluid dynamics, Matematičeskij Sbornik, 47, 271–306, https://hal.archivesouvertes.fr/hal01620642 (last access: 18 December 2023), 1959. a, b, c, d, e, f, g, h
Harris, L., Chen, X., Putman, W., and Zhou, Linjiong andChen, J.H.: A Scientific Description of the GFDL FiniteVolume CubedSphere Dynamical Core, Tech. Rep., National Oceanic and Atmospheric Administration, Office of Oceanic and Atmospheric Research, Geophysical Fluid Dynamics Laboratory (U.S.), https://doi.org/10.25923/6nhs5897, 2021. a
Hertel, O., Berkowicz, R., Christensen, J., and Hoc, O.: Test of two numerical schemes for use in atmospheric transportchemistry models, Atmos. Environ., 27A, 2591–2611, https://doi.org/10.1016/09601686(93)90032T, 1993. a, b, c
Huijnen, V., Williams, J., van Weele, M., van Noije, T., Krol, M., Dentener, F., Segers, A., Houweling, S., Peters, W., de Laat, J., Boersma, F., Bergamaschi, P., van Velthoven, P., Le Sager, P., Eskes, H., Alkemade, F., Scheele, R., Nédélec, P., and Pätz, H.W.: The global chemistry transport model TM5: description and evaluation of the tropospheric chemistry version 3.0, Geosci. Model Dev., 3, 445–473, https://doi.org/10.5194/gmd34452010, 2010. 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
Lachatre, M., Mailler, S., Menut, L., Cholakian, A., Sellitto, P., Siour, G., Guermazi, H., Salerno, G., and Giammanco, S.: Modelling SO_{2} conversion into sulfates in the midtroposphere with a 3D chemistry transport model: the case of Mount Etna's eruption on 12 April 2012, Atmos. Chem. Phys., 22, 13861–13879, https://doi.org/10.5194/acp22138612022, 2022. a, b
Lauritzen, P. H., Skamarock, W. C., Prather, M. J., and Taylor, M. A.: A standard test case suite for twodimensional linear transport on the sphere, Geosci. Model Dev., 5, 887–901, https://doi.org/10.5194/gmd58872012, 2012. a
Lauritzen, P. H., Ullrich, P. A., Jablonowski, C., Bosler, P. A., Calhoun, D., Conley, A. J., Enomoto, T., Dong, L., Dubey, S., Guba, O., Hansen, A. B., Kaas, E., Kent, J., Lamarque, J.F., Prather, M. J., Reinert, D., Shashkin, V. V., Skamarock, W. C., Sørensen, B., Taylor, M. A., and Tolstykh, M. A.: A standard test case suite for twodimensional linear transport on the sphere: results from a collection of stateoftheart schemes, Geosci. Model Dev., 7, 105–145, https://doi.org/10.5194/gmd71052014, 2014. a
Lauritzen, P. H., Conley, A. J., Lamarque, J.F., Vitt, F., and Taylor, M. A.: The terminator “toy” chemistry test: a simple tool to assess errors in transport schemes, Geosci. Model Dev., 8, 1299–1313, https://doi.org/10.5194/gmd812992015, 2015. a, b
LeVeque, R. J.: HighResolution Conservative Algorithms for Advection in Incompressible Flow, SIAM J. Numer. Anal., 33, 627–665, https://doi.org/10.1137/0733033, 1996. a, b, c, d, e, f
Mailler, S.: AdvBench (v1.0.0), Zenodo [code], https://doi.org/10.5281/zenodo.7937121, 2023. a
Mailler, S. and Pennel, R.: ToyCTM v1.0.1, Zenodo [code], https://doi.org/10.5281/zenodo.10018706, 2023. a, b
Mailler, S., Menut, L., di Sarra, A. G., Becagli, S., Di Iorio, T., Bessagnet, B., Briant, R., Formenti, P., Doussin, J.F., GómezAmo, J. L., Mallet, M., Rea, G., Siour, G., Sferlazzo, D. M., Traversi, R., Udisti, R., and Turquety, S.: On the radiative impact of aerosols on photolysis rates: comparison of simulations and observations in the Lampedusa island during the ChArMEx/ADRIMED campaign, Atmos. Chem. Phys., 16, 1219–1244, https://doi.org/10.5194/acp1612192016, 2016. a
Mailler, S., Pennel, R., Menut, L., and Lachâtre, M.: Using the Després and Lagoutière (1999) antidiffusive transport scheme: a promising and novel method against excessive vertical diffusion in chemistrytransport models, Geosci. Model Dev., 14, 2221–2233, https://doi.org/10.5194/gmd1422212021, 2021. a
Mailler, S., Pennel, R., Menut, L., and Cholakian, A.: Scripts to reproduce the results of gmd202378. In Geoscientific Model Development (v1.0.1), Zenodo [code], https://doi.org/10.5281/zenodo.10018761, 2023. a
Martin, R. V., Eastham, S. D., Bindle, L., Lundgren, E. W., Clune, T. L., Keller, C. A., Downs, W., Zhang, D., Lucchesi, R. A., Sulprizio, M. P., Yantosca, R. M., Li, Y., Estrada, L., Putman, W. M., Auer, B. M., Trayanov, A. L., Pawson, S., and Jacob, D. J.: Improved advection, resolution, performance, and community access in the new generation (version 13) of the highperformance GEOSChem global atmospheric chemistry model (GCHP), Geosci. Model Dev., 15, 8731–8748, https://doi.org/10.5194/gmd1587312022, 2022. a
Menut, L., Bessagnet, B., Briant, R., Cholakian, A., Couvidat, F., Mailler, S., Pennel, R., Siour, G., Tuccella, P., Turquety, S., and Valari, M.: The CHIMERE v2020r1 online chemistrytransport model, Geosci. Model Dev., 14, 6781–6811, https://doi.org/10.5194/gmd1467812021, 2021. a
Nair, R. D. and Lauritzen, P. H.: A class of deformational flow test cases for linear transport problems on the sphere, J. Comput. Phys., 229, 8868–8887, https://doi.org/10.1016/j.jcp.2010.08.014, 2010. a
Nakamura, N.: Twodimensional mixing, edge formation, and permeability diagnosed in areacoordinate, J. Atmos. Sci., 53, 1524–1537, https://doi.org/10.1175/15200469(1996)053<1524:TDMEFA>2.0.CO;2, 1996. a, b, c
O'Connor, F. M., Johnson, C. E., Morgenstern, O., Abraham, N. L., Braesicke, P., Dalvi, M., Folberth, G. A., Sanderson, M. G., Telford, P. J., Voulgarakis, A., Young, P. J., Zeng, G., Collins, W. J., and Pyle, J. A.: Evaluation of the new UKCA climatecomposition model – Part 2: The Troposphere, Geosci. Model Dev., 7, 41–91, https://doi.org/10.5194/gmd7412014, 2014. a
Pozzoli, L., Bey, I., Rast, S., Schultz, M. G., Stier, P., and Feichter, J.: Trace gas and aerosol interactions in the fully coupled model of aerosolchemistryclimate ECHAM5HAMMOZ: 1. Model description and insights from the spring 2001 TRACEP experiment, J. Geophys. Res.Atmos., 113, D07308, https://doi.org/10.1029/2007JD009007, 2008. a
Prather, M. J.: Numerical advection by conservation of secondorder moments, J. Geophys. Res.Atmos., 91, 6671–6681, https://doi.org/10.1029/JD091iD06p06671, 1986. a
Putman, W. M. and Lin, S.J.: Finitevolume transport on various cubedsphere grids, J. Comput. Phys., 227, 55–78, https://doi.org/10.1016/j.jcp.2007.07.022, 2007. a
Seinfeld, J. H. and Pandis, S. N.: Atmospheric chemistry and physics: From air pollution to climate change, WileyInterscience, ISBN 0471178152, 1997. a
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
Timmermans, R., van Pinxteren, D., Kranenburg, R., Hendriks, C., Fomba, K., Herrmann, H., and Schaap, M.: Evaluation of modelled LOTOSEUROS with observational based PM_{10} source attribution, Atmospheric Environment: X, 14, 100173, https://doi.org/10.1016/j.aeaoa.2022.100173, 2022. 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
Walcek, C. J.: Minor flux adjustment near mixing ratio extremes for simplified yet highly accurate monotonic calculation of tracer advection, J. Geophys. Res., 105, 9335–9348, https://doi.org/10.1029/1999JD901142, 2000. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s
Waruszewski, M., Kühnlein, C., Pawlowska, H., and Smolarkiewicz, P. K.: MPDATA: Thirdorder accuracy for variable flows, J. Comput. Phys., 359, 361–379, https://doi.org/10.1016/j.jcp.2018.01.005, 2018. a
Zhao, S., Russell, M. G., Hakami, A., Capps, S. L., Turner, M. D., Henze, D. K., Percell, P. B., Resler, J., Shen, H., Russell, A. G., Nenes, A., Pappin, A. J., Napelenok, S. L., Bash, J. O., Fahey, K. M., Carmichael, G. R., Stanier, C. O., and Chai, T.: A multiphase CMAQ version 5.0 adjoint, Geosci. Model Dev., 13, 2925–2944, https://doi.org/10.5194/gmd1329252020, 2020. a