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

This study presents a novel method to estimate the performance of advection schemes in numerical experiments along with a semi-realistic, non-linear, 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 chemistry-transport modellers an alternative, high-performance scheme designed for Cartesian-grid Eulerian chemistry-transport 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.

1 Introduction

Chemistry-transport 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 chemistry-transport models are the flux-based 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 chemistry-transport models. For example, GEOS-Chem 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 chemistry-transport models including CCATT-BRAMS (Freitas et al., 2012) and LOTOS-EUROS (Timmermans et al., 2022). Therefore, while not the only ones used in chemistry-transport models, the schemes we study here are the base of the numerical resolution of advection in some of the most common chemistry-transport models. Other popular schemes include versions of the Bott (1989) scheme, which is less diffusive than PPM but has the inconvenience of being non-monotonic 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 higher-order 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 chemistry-transport 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 chemistry-transport 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 above-cited 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 non-linear chemical system, they have tested chemistry-advection combinations, and the errors they generate as diagnostics of mass-conservation 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 non-linear chemical processes in the troposphere (in their case, the oxidation pathways of SO₂ in a mid-tropospheric 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 chemistry-transport 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 ECHAM5-HAMMOZ model (Pozzoli et al., 2008), the TM5 model (Huijnen et al., 2010) and the UKCA climate-composition 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 time-stepping 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 Numerical-experiment description

2.1 Chemical mechanism

The chemical mechanism used here (Reactions R1–R12) includes a subset of the main reactions of tropospheric gas-phase 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₂ by HO₂. 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.

$$\begin{array}{}\text{(R1)}& {\displaystyle }& {\displaystyle }{\mathrm{NO}}_{\mathrm{2}}+h\mathit{\nu }\to \mathrm{NO}+\mathrm{O}\text{(R2)}& {\displaystyle }& {\displaystyle }\mathrm{O}+{\mathrm{O}}_{\mathrm{2}}+M\to {\mathrm{O}}_{\mathrm{3}}+M\text{(R3)}& {\displaystyle }& {\displaystyle }\mathrm{NO}+{\mathrm{O}}_{\mathrm{3}}\to {\mathrm{NO}}_{\mathrm{2}}+{\mathrm{O}}_{\mathrm{2}}\text{(R4)}& {\displaystyle }& {\displaystyle }{\mathrm{O}}_{\mathrm{3}}+h\mathit{\nu }\to \mathrm{O}{(}^{\mathrm{1}}\mathrm{D})+{\mathrm{O}}_{\mathrm{2}}\text{(R5)}& {\displaystyle }& {\displaystyle }\mathrm{O}{(}^{\mathrm{1}}\mathrm{D})+{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\to \mathrm{OH}+\mathrm{OH}\text{(R6)}& {\displaystyle }& {\displaystyle }\mathrm{O}{(}^{\mathrm{1}}\mathrm{D})+{\mathrm{N}}_{\mathrm{2}}\to \mathrm{O}+{\mathrm{N}}_{\mathrm{2}}\text{(R7)}& {\displaystyle }& {\displaystyle }\mathrm{O}{(}^{\mathrm{1}}\mathrm{D})+{\mathrm{O}}_{\mathrm{2}}\to \mathrm{O}+{\mathrm{O}}_{\mathrm{2}}\text{(R8)}& {\displaystyle }& {\displaystyle }\mathrm{CO}+\mathrm{OH}\to {\mathrm{CO}}_{\mathrm{2}}+{\mathrm{HO}}_{\mathrm{2}}\text{(R9)}& {\displaystyle }& {\displaystyle }\mathrm{NO}+{\mathrm{HO}}_{\mathrm{2}}\to {\mathrm{NO}}_{\mathrm{2}}+\mathrm{OH}\text{(R10)}& {\displaystyle }& {\displaystyle }{\mathrm{HO}}_{\mathrm{2}}+{\mathrm{HO}}_{\mathrm{2}}\to {\mathrm{H}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{2}}+{\mathrm{O}}_{\mathrm{2}}\text{(R11)}& {\displaystyle }& {\displaystyle }\mathrm{OH}+{\mathrm{HO}}_{\mathrm{2}}\to {\mathrm{H}}_{\mathrm{2}}\mathrm{O}+{\mathrm{O}}_{\mathrm{2}}\text{(R12)}& {\displaystyle }& {\displaystyle }{\mathrm{NO}}_{\mathrm{2}}+\mathrm{OH}\to {\mathrm{HNO}}_{\mathrm{3}}\end{array}$$

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}}{\mathrm{s}}^{-\mathrm{1}}$

• $j_{\mathrm{R}\mathrm{4}}=\mathrm{2.5}\times {\mathrm{10}}^{-\mathrm{5}}{\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 chemistry-transport 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 pseudo-reaction rate constant that depends on the concentration of air molecules and of water vapour molecules and on applying the quasi-steady-state approximation to O(¹D). For this study, we have chosen to treat O(¹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(¹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;

• non-linear 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⁵ 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 regional-scale chemistry-transport modelling. Domain 𝒟 is therefore discretized into n² 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):

$$\begin{array}{}\text{(1)}& \begin{array}{rl}u& ={\displaystyle \frac{L}{T}}{\sin }^{\mathrm{2}}\left({\displaystyle \frac{\mathit{\pi }x}{L}}\right)\sin \left(\mathrm{2}\mathit{\pi }y\right)g\left(t\right),\\ & v=-{\displaystyle \frac{L}{T}}{\sin }^{\mathrm{2}}\left({\displaystyle \frac{\mathit{\pi }y}{L}}\right)\sin \left(\mathrm{2}\mathit{\pi }x\right)g\left(t\right),\end{array}\end{array}$$

with

$$\begin{array}{}\text{(2)}& g\left(t\right)=\cos \left({\displaystyle \frac{\mathit{\pi }t}{T}}\right).\end{array}$$

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 non-dimensional scales for time and space, we set a dimensional scale length L=10⁵ 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 time-dependent streamfunction for this flow is

$$\begin{array}{}\text{(3)}& \mathit{\psi }\left(x,y,t\right)=-{\displaystyle \frac{L^{\mathrm{2}}}{\mathit{\pi }T}}{\sin }^{\mathrm{2}}\left({\displaystyle \frac{\mathit{\pi }x}{L}}\right){\sin }^{\mathrm{2}}\left({\displaystyle \frac{\mathit{\pi }y}{L}}\right)g\left(t\right).\end{array}$$

Figure 1 Streamlines (black contours), wind vectors and wind module in m s⁻¹ (colour shades) at t=0 for the swirling deformational wind field defined in Eq. (1) with L=10⁵ m and T=86 400 s.

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:

$$\begin{array}{}\text{(4)}& {\displaystyle }\mathit{\phi }\left(x,y\right)& {\displaystyle }={\sin }^{\mathrm{2}}{\displaystyle \frac{\mathrm{2}\mathit{\pi }x}{L}}{\sin }^{\mathrm{2}}{\displaystyle \frac{\mathrm{2}\mathit{\pi }y}{L}}\mathrm{if}x<{\displaystyle \frac{L}{\mathrm{2}}}\text{and}x<{\displaystyle \frac{L}{\mathrm{2}}},\text{(5)}& {\displaystyle }\mathit{\phi }\left(x,y\right)& {\displaystyle }=\mathrm{0}\text{otherwise}.\end{array}$$

The initial conditions are defined as follows (in terms of the mixing ratio):

• ${\mathit{\alpha }}_{\mathrm{TRC}}=\mathrm{100}\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}\mathrm{ppb}\times \mathit{\phi }\left(x,y\right)$

• ${\mathit{\alpha }}_{{\mathrm{NO}}_{\mathrm{2}}}=\mathrm{10}\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⁻¹. The initial concentrations of the other active species (O, O(¹D), OH, HO₂, CO₂, H₂O₂, HNO₃) are initialized to zero. They will be produced by Reactions (R1), (R4), (R5), (R8) and (R12).

Figure 2 Initial mixing ratio of TRC, proportional to φ(x,y) (defined in Eqs. 4–5), discretized on domain 𝒟 with n=25 subintervals ($\mathit{\delta }x=\mathit{\delta }y=\mathrm{4000}$ m). Note that, at t=0, α_NO=α_TRC and ${\mathit{\alpha }}_{{\mathrm{NO}}_{\mathrm{2}}}=\frac{{\mathit{\alpha }}_{\mathrm{TRC}}}{\mathrm{10}}$ have the same spatial distribution.

Finally, another species of inert tracer, TRCb, is introduced so that

$$\begin{array}{}\text{(6)}& {\mathit{\alpha }}_{\mathrm{TRCb}}+{\mathit{\alpha }}_{\mathrm{NO}}+{\mathit{\alpha }}_{{\mathrm{NO}}_{\mathrm{2}}}=\mathrm{110}\mathrm{ppb}.\end{array}$$

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 Numerical methods

3.1 Advection schemes

3.1.1 Existing advection schemes

The following existing advection schemes have been tested in the study:

1. Godunov (1959)

2. Van Leer (1977)

3. Walcek (2000)

4. PPM (Colella and Woodward, 1984).

These schemes are flux-based, upwind-biased, semi-Lagrangian schemes based on polynomial reconstructions of the average concentrations. These polynomial reconstructions are piecewise-constant for Godunov (1959), piecewise-linear 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 upwind-biased 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 mixing-ratio 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{̃}}{\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{̃}}{\mathit{\alpha }}}_{i+\frac{\mathrm{1}}{\mathrm{2}}}$ following the piecewise parabolic method procedure described in Colella and Woodward (1984).

• Otherwise we estimate the Walcek-adjusted flux as follows:
  $s=\mathrm{sign}\left({\mathit{\alpha }}_{i+\mathrm{1}}-{\mathit{\alpha }}_i\right)\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 }=\max \left(\mathrm{1.5},\mathrm{1.2}+\mathrm{0.6}\mathit{\nu }\right)$,
  ${\stackrel{\mathrm{̃}}{\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.

Figure 3 Reconstruction of a Gaussian mixing-ratio profile by (a) the Van Leer (1977) scheme, (b) the Walcek (2000) scheme, (c) the PPM scheme and (d) the PPM + W scheme. The x axis is a non-dimensional space coordinate. The reconstruction has been performed for a Courant number ν=0.4. The reconstructed fields are presented with alternating red and blue colours to enhance the discontinuities between neighbouring cells.

Figure 4 (a) Same as Fig. 3c but zoomed in to within the vicinity of the mixing-ratio maximum and (b) same as Fig. 3d but zoomed in to within the vicinity of the mixing-ratio maximum.

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⁵ cells has been performed over 520 time steps, corresponding to 1.04×10⁸ calls to the reconstruction routine, plus the update of the mixing-ratio 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).

Table 1 Mean calculation time per cell and per time step for the five advection schemes retained for the present study. The calculation has been performed in Fortran, a programming language frequently used for operational chemistry-transport models, on a laptop with an Intel Core i7-1165G7 CPU.

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 ∥⋅∥₁ (Fig. 5a) and ∥⋅∥₂ (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 ∥⋅∥₁ convergence rate and around 2 for ∥⋅∥₂ (Fig. 2).

Figure 5 (a) The ∥⋅∥₁ error as a function of the number of n_x for the five advection schemes used in the present study; (b) same as (a) for the ∥⋅∥₂ error.

Table 2 Convergence rate of the five advection schemes used in the present study for ∥⋅∥₁ and ∥⋅∥₂.

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

$$\begin{array}{}\text{(7)}& \mathit{c}\left(t+\mathit{\delta }{t}_{\mathrm{chem}}\right)={\displaystyle \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)}}.\end{array}$$

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 non-linear, fixed-point 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 chemistry-transport 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 best-possible 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⁻¹ (Fig. 1) and $\mathit{\delta }x=\mathrm{4}\times {\mathrm{10}}^{\mathrm{3}}$ m, this yields a maximal Courant number ν_max≃0.8. The time-stepping strategy follows Strang-style time stepping (Strang, 1968), with the steps as follows:

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

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

3. Integrate meridional advection over δt_adv.

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

5. 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 above-described 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)

Table 3 Summary of the main characteristics of the simulations that have been performed.

3.4 Conservation properties

It is worth noting that, by construction, flux-based advection integration is mass-conservative 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₂) and H can be expected to be of the order of ε, since the Euler backward iterative scheme (Eq. 7) is in principle mass-conservative, 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 flux-based 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.

Table 4 Mass-conservation diagnostic for C, active N, H and TRC in the PPM + W and base simulations.

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 mixing-ratio extrema.

4 The signature function

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 ∥⋅∥₁, introducing E₁ as the normalized $\Vert \cdot {\Vert }_{\mathrm{1}}$ error in the mixing ratio:

$$\begin{array}{}\text{(8)}& E_{\mathrm{1}}={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^N{\mathit{\rho }}_i{V}_i\left|{\mathit{\alpha }}_i^{\mathrm{t}}-{\mathit{\alpha }}_i^{\mathrm{0}}\right|}{{\sum }_{i=\mathrm{1}}^N{\mathit{\rho }}_i{V}_i{\mathit{\alpha }}_i^{\mathrm{0}}}},\end{array}$$

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

$$\begin{array}{}\text{(9)}& E_{\mathrm{1}}={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^N\left|{\mathit{\alpha }}_i^{\mathrm{t}}-{\mathit{\alpha }}_i^{\mathrm{0}}\right|}{{\sum }_{i=\mathrm{1}}^N{\mathit{\alpha }}_i^{\mathrm{0}}}}.\end{array}$$

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 area-coordinate formulation used by Nakamura (1996).

Let us imagine a fluid with density $\mathit{\rho }\left(x;y;z;t\right)$ in a three-dimensional 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

$$\begin{array}{}\text{(10)}& {\displaystyle \frac{\partial \mathit{\rho }}{\partial t}}+\mathrm{\nabla }\left(\mathit{\rho }\mathit{u}\right)=\mathrm{0},\text{(11)}& {\displaystyle \frac{\partial \mathit{\alpha }}{\partial t}}+\mathit{u}\cdot \mathrm{\nabla }\mathit{\alpha }=\mathrm{0}.\end{array}$$

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

$$\begin{array}{}\text{(12)}& {\mathcal{S}}^{\mathrm{t}}\left(X\right)={\displaystyle \frac{{\int }_{\mathcal{D}}H\left(X-\mathit{\alpha }\left(x;y;z;t\right)\right)\mathit{\rho }\mathrm{d}\mathcal{V}}{{\int }_{\mathcal{D}}\mathit{\rho }\mathrm{d}\mathcal{V}}},\end{array}$$

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 mass-weighted 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 area-coordinate 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, 𝒮⁰(X) is the (normalized) mass of the fluid parcel 𝒟⁰(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 𝒟⁰(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)=𝒮⁰(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 non-divergent flows.

Therefore, since we know that, for the exact solution, 𝒮^t=𝒮⁰, the departure of the numerical evaluation of function 𝒮^t from 𝒮⁰ 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:

$$\begin{array}{}\text{(13)}& {\mathcal{S}}^{\mathrm{t}}\left(X\right)\simeq {\displaystyle \frac{{\sum }_{i=\mathrm{1}}^{N}H\left(X-{\mathit{\alpha }}_i\right){\mathit{\rho }}_i{V}_i}{{\sum }_{i=\mathrm{1}}^N{\mathit{\rho }}_i{V}_i}}.\end{array}$$

Generally, the numerical evaluation of 𝒮^t in an Eulerian model will differ from 𝒮⁰: 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 non-divergent and the carrier fluid mass ρ_iV_i is the same in all model cells, the norm-1 difference between 𝒮^t and 𝒮⁰ can be calculated as

$$\begin{array}{}\text{(14)}& \underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\left|{\mathcal{S}}^{\mathrm{t}}\left(X\right)-{\mathcal{S}}^{\mathrm{0}}\left(X\right)\right|\mathrm{d}X={\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^N\left|{\stackrel{\mathrm{̃}}{{\mathit{\alpha }}^{\mathrm{t}}}}_i-{\stackrel{\mathrm{̃}}{{\mathit{\alpha }}^{\mathrm{0}}}}_i\right|,\end{array}$$

where ${\stackrel{\mathrm{̃}}{{\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 norm-1 difference ℰ₁ between 𝒮^t and 𝒮⁰ as

$$\begin{array}{}\text{(15)}& {\mathcal{E}}_{\mathrm{1}}={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^N\left|{\stackrel{\mathrm{̃}}{{\mathit{\alpha }}^{\mathrm{t}}}}_i-{\stackrel{\mathrm{̃}}{{\mathit{\alpha }}^{\mathrm{0}}}}_i\right|}{{\sum }_{i=\mathrm{1}}^N{\stackrel{\mathrm{̃}}{{\mathit{\alpha }}^{\mathrm{0}}}}_i}}.\end{array}$$

Figure 6 shows an example of how comparing the 𝒮^t signature function with 𝒮⁰ 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 𝒮⁰ for the Godunov and Van Leer simulations, respectively. This graphical representation permits giving an intuitive meaning to the normalized signature error ℰ₁ as the total area between the representative curves of 𝒮^t and 𝒮⁰ (shaded in Fig. 6c–d), divided by the total area left of the representative curve of 𝒮⁰ (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 𝒮⁰ and 𝒮^t is smaller in the Van Leer simulation (Fig. 6d, with ℰ₁=0.218) than in the Godunov simulation (Fig. 6c, with ℰ₁=0.571).

Figure 6 (a) The TRC mixing ratio as simulated in the Godunov simulation at time $\frac{T}{\mathrm{2}}$; (b) same as (a) for the Van Leer simulation. (c) The signature function for the TRC mixing ratio in the Godunov simulation compared to the base simulation (𝒮⁰, green line); (d) same as (c) for the Van Leer simulation. The signature error ℰ₁ is equal to the ratio of the shaded area between the representative curves of 𝒮⁰ and 𝒮^t to the hatched area left of 𝒮⁰.

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{‾}}{\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{‾}}{\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 non-linear chemical system such as Reactions (R1)–(R12), there is no easy access to the exact solution of the system even without advection: ${\stackrel{\mathrm{‾}}{\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{‾}}{\mathcal{S}}}^{\mathrm{t}}$ and the objective calculation of ℰ₁ clearly show that the PPM + W simulation agrees better with the base simulation in terms of the signature function (ℰ₁=0.0622), followed by the PPM simulation (ℰ₁=0.0679) and the Van Leer simulation (ℰ₁=0.0787). Apart from this quantitative agreement, some qualitative and local conclusions can be drawn from the graphical comparisons between 𝒮^t and ${\stackrel{\mathrm{‾}}{\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.

Figure 7 Ozone mixing ratio at $\frac{T}{\mathrm{2}}$ as simulated in (a) the Van Leer simulation, (b) the PPM simulation, (c) the PPM + W simulation and (d) the base simulation.

Figure 8 The signature function for the O₃ mixing ratio at $\frac{T}{\mathrm{2}}$, ${\mathcal{S}}^{\frac{T}{\mathrm{2}}}$, as simulated in (a) the Van Leer simulation; (b) the PPM simulation and (c) the PPM + W simulation. As in Fig. 6, the signature error ℰ₁ is equal to the ratio of the shaded area between the representative curves of 𝒮⁰ and 𝒮^t to the hatched area left of 𝒮⁰.

5 Results and discussion

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₃, 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 non-dimensional 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.

Figure 9 Time series for all the simulations for the (a) minimum and maximum of the O₃ mixing ratio, (b) maximum TRC mixing ratio, and (c) half volume of the TRC plume (relative to the volume of the entire domain).

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 ℰ₁ 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 ℰ₁ for TRC, O₃, NO and NO₂ for all the simulations, with the smallest ℰ₁ indicating the best-performing 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 ℰ₁ 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 ℰ₁ throughout the simulation shows that, for TRC, NO and NO₂, the best-performing 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.

Figure 10 Time series for ℰ₁ ($\Vert \cdot {\Vert }_{\mathrm{1}}$ normalized signature error) for (a) TRC, (b) O₃, (c) NO and (d) NO₂.

Table 5 Normalized $\Vert \cdot {\Vert }_{\mathrm{1}}$ error E₁ (Eq. 9) and normalized $\Vert \cdot {\Vert }_{\mathrm{1}}$ signature error ℰ₁ (Eq. 15) at the end of the simulations for O₃, NO, NO₂ and TRC, compared to the base simulation with no advection. In each column, the lowest error value is in bold font the second-lowest is underlined.

While we have shown so far that the analysis of ℰ₁ permits the drawing of clearer conclusions regarding the performance of the various simulations, in the present case we can confirm these results by comparing ℰ₁ to a more classical metric, the normalized $\Vert \cdot {\Vert }_{\mathrm{1}}$ error E₁ (Eq. 9) of the simulations with advection relative to the base simulation without advection. Unlike ℰ₁, E₁ 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₁ and ℰ₁ in Table 5.

First, in all cases we always have ℰ₁<E₁. Qualitatively, this can be interpreted as ℰ₁ 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 ℰ₁ is the same for all variables as with E₁: for TRC, NO and NO₂, PPM + W performs best, followed in this order by Walcek, PPM, Van Leer and Godunov, but for O₃, PPM performs better than Walcek. In our study case, analysing ℰ₁ 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₁ (which requires access to the exact solution). Even for simulations with comparable performance to the Walcek and PPM simulations, the signature error ℰ₁ 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 ℰ₁ are the same and correspond to E₁. This gives confidence in the usability of the signature error ℰ₁ 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 ℰ₁ 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₂, 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₃, having an initially uniform distribution, errors due to advection appear only when sufficient heterogeneity is introduced into the O₃ 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₃ 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 higher-resolution 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₂. However, PPM + W clearly performs better for O₃ 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).

Table 6 Normalized $\Vert \cdot {\Vert }_{\mathrm{1}}$ error E₁ and normalized $\Vert \cdot {\Vert }_{\mathrm{1}}$ signature error ℰ₁ in the 1 km simulations at the end of the simulations for O₃, NO, NO₂ and TRC, compared to the base simulation at 1 km resolution with no advection. In each column, the lowest error value is in bold font and the second-lowest is underlined.

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, second-order, monotonic schemes for the advection equation are all non-linear. In particular, among all the schemes we study here, all are non-linear due to their specific treatment of the maxima to preserve monotonicity. However, these non-linear 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 higher-order and monotonic, but therefore non-linear, 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 non-linearity of advection schemes, which is an undesirable property since the advection equation itself is linear. Figure 11 shows that non-linearities of up to 1 %–2 % appear for all schemes in the representation of CO+CO₂ 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 non-linearities reach up to 10 % of the expected value and are strongest in the PPM + W scheme. These results suggest that non-conservation of sum species in the vicinity of mixing-ratio extrema, which is observed in all second-order 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 non-linear behaviour in the vicinity of the extrema, PPM + W preserves the values of extrema better than PPM (Fig. 9b).

Figure 11 Minimum and maximum values in the 4 km simulations (presented in the paper) for ${\mathit{\alpha }}_{\mathrm{CO}}+{\mathit{\alpha }}_{{\mathrm{CO}}_{\mathrm{2}}}$ (panel a) and ${\mathit{\alpha }}_{\mathrm{NO}}+{\mathit{\alpha }}_{{\mathrm{NO}}_{\mathrm{2}}}+{\mathit{\alpha }}_{{\mathrm{HNO}}_{\mathrm{3}}}+{\mathit{\alpha }}_{\mathrm{TRCb}}$ (panel b).

6 Conclusions

We have introduced the signature function 𝒮^t as a sort of mass-weighted 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 above-cited, 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 non-conservation of 𝒮^t gives a proxy for model error. Based on this idea, we propose the normalized $\Vert \cdot {\Vert }_{\mathrm{1}}$ signature error ℰ₁ (Eq. 15) as a measure of model error. Graphically, ℰ₁ 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 ℰ₁ 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 ℰ₁ 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 short-lived species such as O(¹D) or O (Reactions R1–R12). With an initial peak of NO and NO₂ concentrated over a small area (Fig. 2) and initially uniform concentrations of CO and O₃, this case generates a sharp ozone minimum colocalized with the NO_x peak, a large area with background O₃ 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₁ of 16.4 % in O₃ (19.3 % for PPM), 20.7 % in TRC (24.3 % for Walcek), 28.3 % in NO₂ (31.2 % for Walcek) and 37.2 % in NO (42.5 % for Walcek). Table 5 also shows the ℰ₁ signature error, which is always smaller than E₁. Interestingly, examining the ℰ₁ 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₁. This shows that, even without access to an exact solution, the ℰ₁ 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₁. This being shown, Fig. 10 permits visualization of the evolution of error along the course of the simulation, while E₁ 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 chemistry-transport models, generalizing this concept to more dimensions (by studying the mass-weighted 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 ℰ₁ 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 real-world geophysical models is that it relies on mixing-ratio 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₂), PPM + W is even outperformed by PPM in terms of accuracy. However, even in this high-resolution 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 non-linearity 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 chemistry-transport 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 chemistry-transport models, in an effort to reduce numerical diffusion, which is important in particular in the presence of non-linear chemistry as discussed in Lachatre et al. (2022).