A novel Eulerian model based on central moments to simulate age and reactivity continua interacting with mixing processes

In geoscientific models, simulating the properties associated with particles in a continuum can serve many scientific purposes, and this has commonly been addressed using Lagrangian models. As an alternative approach, we present an Eulerian method here: diffusion–advection–reaction type partial differential equations are derived for centralized moments, which can describe the distribution of properties associated with chemicals in reaction–transport models. When the property is age, the equations for centralized moments (unlike non-central moments) do not require terms to account for aging, making this method suitable for modeling age tracers. The properties described by the distributions may also represent kinetic variables affecting reaction rates. In practical applications, continuous distributions of ages and reactivities are resolved to simulate organic matter mineralization in surficial sediments, where macrofaunal and physical mixing processes typically dominate transport. In test simulations, mixing emerged as the predominant factor shaping reactivity and age distributions. Furthermore, the applications showcase the method's aptitude for modeling continua in mixed environments while also highlighting practical considerations and challenges.

1 Introduction

The partial differential equation (PDE) to describe chemical diffusion (Fick, 1855) is mathematically equivalent to Fourier's heat conduction equation (Fourier, 1822), which has become ubiquitous in science to describe random transport processes (Narasimhan, 1999). When materials are transported, tracking other associated properties besides the concentration can be desirable. For example, the age (or transit time) of fluids and chemicals has often been simulated. Diffusive mixing will lead to a local spreading in ages. Resolving these distributions is commonly achieved through Lagrangian approaches, which aim to simulate sufficiently large numbers of individual particles to describe the evolution of the properties and their statistical distributions. Alternatively, Eulerian approaches utilizing PDEs offer the advantage of analytical evaluation and are computationally less expensive. Analytical solutions for age distributions and boundary condition problems in particular can be found in Delhez and Deleersnijder (2002) and Kuderer (2022). Deleersnijder et al. (2001) and Delhez and Deleersnijder (2002) derived Eulerian PDEs to simulate the effect of diffusion on the mean and higher raw moments and considered the effect of radioactive decay on age distributions. In this study, we derive Eulerian PDEs for centralized moments. These are more readily intuitively understood than raw moments and are not affected by aging, making them ideal for modeling time tracers.

Beyond modeling passive tracers, we intend to test moment-based PDEs in more complex applications whereby chemical reactions depend on and affect distributions. Modeling the effect of “aging” on the apparent organic matter reactivity (Middelburg, 1989, 2019) provides an interesting practical case study. Bulk organic matter in sediments and soils contains materials with varying reactivities (e.g., De Leeuw and Largeau, 1993), and the bulk degradation rate depends on the entire matrix. As more reactive components disappear first, the remaining organic matter becomes more refractory (Zonneveld et al., 2010). Organic molecules also undergo transformations, which generally lower the reactivity (Burdige, 2007). The overall decreasing reactivity over time is contained in the concept of aging. In multi-G models, separate state variables represent discrete classes of varying reactivities (Jørgensen, 1978; Westrich and Berner, 1984). The disadvantage of this approach is that reactivity classes and their distribution in deposited organic matter are somewhat arbitrarily chosen (Jørgensen, 1978), resulting in parameterizations that are difficult to compare between studies. Also, no more than three classes are usually defined, which cannot represent more gradual changes in reactivity. The reactivity of organic matter may also be described as a continuum for which various distribution functions have been proposed (e.g., Boudreau and Ruddick, 1991; Vähätalo et al., 2010; Xu et al., 2023). The gamma distribution is most commonly used (Arndt et al., 2013; Freitas et al., 2021), in part because it allows an analytical solution for the evolution of the continuum over time (Boudreau and Ruddick, 1991). It can be easily implemented in sediment models by replacing time with sediment depth based on the assumption of a constant burial velocity or a reconstruction of the deposition history. However, the space-for-time substitution only accounts for the burial of particulate organic matter and ignores mixing processes.

Animals and plants continuously cause disturbances in sediments, which is referred to as bioturbation in the literature (Meysman et al., 2006). Bioturbation typically dominates the transport of solids in sediments up to a depth of ∼ 10 cm (Tromp et al., 1995; Middelburg et al., 1997; Boudreau, 1994). In reaction–transport models, this process is most commonly implemented as Fickian transport (Goldberg and Koide, 1962; Guinasso and Schink, 1975; Meysman et al., 2005), i.e., as chemical diffusion, but with a diffusivity decreasing over depth. Mixing of particulate organic matter can also be the result of other natural processes or anthropogenic activities, such as trawl fishing (e.g., De Borger et al., 2021). The inability of reactivity continuum models to account for mixing either caused an error or limited their application to environments without significant mixing processes (Freitas et al., 2021). Previously, Lagrangian methods have been developed to simulate organic matter mineralization in turbated sediments (Meile and Van Cappellen, 2005; Kuderer, 2022), but these have only included reaction networks with few chemicals. The development of an alternative Eulerian approach, compatible with classical early diagenetic reaction–transport models (Wang and Van Cappellen, 1996; Boudreau, 1997), may be needed to encourage wider usage of modeled reactivity continua in turbated environments.

Here, we first derive expressions that describe the effect of diffusion on central moments in partial differential equations (Sect. 2). Next, we develop an approach to derive additional reaction terms that may depend on the distribution (Sect. 3), expanding beyond previous studies that used raw moments to simulate transit time distributions and only considered radioactive decay. Finally, the central-moment-based method is tested for age distributions and organic matter mineralization in bioturbated sediment (Sect. 4) to explore whether a moments-based approach can effectively describe age or reactivity continua in reaction–transport models.

2 Derivation of partial differential equations for diffusion

Diffusion is a process that mixes distributions of properties associated with moving particles. In the derivation, we will assume that the property of interest is age, even though it could be any other scalar property that does not affect transport. First, we derive equations for chemical diffusion (see Sect. 2.1) and the effect of diffusion on mean age (see Sect. 2.2) to illustrate the method based on microscopic diffusion. We then derive partial differential equations for higher centralized moments (see Sect. 2.3).

2.1 Microscopic derivation for concentration

Following Crank (1956), microscopic diffusion can be represented as random jumps back and forth. Consider three locations, left (L), center (C), and right (R), aligned on a line and separated by the jumping distance of particles δ_x. The change in the number of molecules at location “C” is given by

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\mathrm{\Delta }{n}_{\mathrm{C}}}{\mathrm{\Delta }t}}=\mathrm{0.5}{f}_{\mathrm{r}}(n_{\mathrm{R}}-n_{\mathrm{C}})-\mathrm{0.5}{f}_{\mathrm{l}}(n_{\mathrm{C}}-n_{\mathrm{L}}),\end{array}$$

where f_l and f_r are the jumping frequencies in the leftward and rightward directions, respectively, and n_X is the number of particles at location X. Lower-case and upper-case subscripts indicate evaluations at the boundaries and centers of cells, respectively. After dividing by volume V, defining $C=n/V$, and also multiplying the right-hand side by ${\mathit{\delta }}_x^{\mathrm{2}}/{\mathit{\delta }}_x^{\mathrm{2}}$,

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\mathrm{\Delta }{C}_{\mathrm{C}}}{\mathrm{\Delta }t}}={\displaystyle \frac{\mathrm{1}}{{\mathit{\delta }}_x}}\left({\displaystyle \frac{\mathrm{0.5}{f}_{\mathrm{r}}{\mathit{\delta }}_x^{\mathrm{2}}(C_{\mathrm{R}}-C_{\mathrm{C}})}{{\mathit{\delta }}_x}}-{\displaystyle \frac{\mathrm{0.5}{f}_{\mathrm{l}}{\mathit{\delta }}_x^{\mathrm{2}}(C_{\mathrm{C}}-C_{\mathrm{L}})}{{\mathit{\delta }}_x}}\right)\end{array}$$

is obtained. The diffusivity is identified as $D=\mathrm{0.5}f{\mathit{\delta }}_x^{\mathrm{2}}$. One linearization is made, i.e., $\mathrm{\Delta }C=\partial C/\partial x{\mathit{\delta }}_x$, resulting in

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\mathrm{\Delta }{C}_{\mathrm{C}}}{\mathrm{\Delta }t}}={\displaystyle \frac{\mathrm{1}}{{\mathit{\delta }}_x}}\left(D_{\mathrm{r}}{\left.{\displaystyle \frac{\partial C}{\partial x}}\right|}_{\mathrm{r}}-D_{\mathrm{l}}{\left.{\displaystyle \frac{\partial C}{\partial x}}\right|}_{\mathrm{l}}\right).\end{array}$$

Applying the divergence theorem yields

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\partial C}{\partial t}}={\displaystyle \frac{\partial }{\partial x}}\left(D{\displaystyle \frac{\partial C}{\partial x}}\right),\end{array}$$

which is the diffusion equation.

Table 1 Notation.

2.2 Microscopic derivation of the diffusion equation for the mean

The same method is applied to the mean age associated with particles, which is the first raw moment. Let χ_i be the age of a particle and ${\sum }_{i=\mathrm{1}}^n{\mathit{\chi }}_i$ the total age of all particles n in a control volume V, meaning that the mean age of the particles is $\mathit{\mu }={\sum }_i{\mathit{\chi }}_i/n$. As a result, $\mathit{\mu }C={\sum }_i{\mathit{\chi }}_i/V$ is the summed ages of all particles per control volume.

Let j₁ and j₄ be the fluxes that transport particles into the control volume from left and right, respectively. Similarly, let j₂ and j₃ be the fluxes that remove matter in the rightward and leftward direction, respectively. Substituting the summed total age of particles for the total number of particles in Eq. (1) yields

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{\mathrm{\Delta }\left(C_{\mathrm{C}}{\mathit{\mu }}_{\mathrm{C}}\right)}{\mathrm{\Delta }t}}={\displaystyle \frac{\mathrm{1}}{V}}\left[j_{\mathrm{4}}{\mathit{\mu }}_{j_{\mathrm{4}}}-j_{\mathrm{3}}{\mathit{\mu }}_{j_{\mathrm{3}}}-(j_{\mathrm{2}}{\mathit{\mu }}_{j_{\mathrm{2}}}-j_{\mathrm{1}}{\mathit{\mu }}_{j_{\mathrm{1}}})\right],\end{array}$$

whereby μ_j is the mean age of the jumping particles, and the fluxes j_k have dimensions using the number of particles over time. Note that this section only considers changes to the local mean age caused by diffusive transport. In Sect. 3.2 derivatives to account for the effect of aging will be discussed. From Eq. (1) it follows that j₁=0.5f_ln_L, j₂=0.5f_ln_C, j₃=0.5f_rn_C, and j₄=0.5f_rn_R. When it is assumed that the random jumps are not affected by age, the mean age of a larger number of jumping particles will approach the mean age at the source location X, meaning that $\langle {\mathit{\mu }}_{j_k}\rangle ={\mathit{\mu }}_X$. Making this substitution and repeating the steps that were taken for the derivation of chemical diffusion yields

$$\begin{array}{}\text{(6)}& \begin{array}{rl}& 〈{\displaystyle \frac{\mathrm{\Delta }\left(C_{\mathrm{C}}{\mathit{\mu }}_{\mathrm{C}}\right)}{\mathrm{\Delta }t}}〉={\displaystyle \frac{\mathrm{1}}{{\mathit{\delta }}_x}}\\ & \left({\displaystyle \frac{D_{\mathrm{r}}(C_{\mathrm{R}}{\mathit{\mu }}_{\mathrm{R}}-C_{\mathrm{C}}{\mathit{\mu }}_{\mathrm{C}})}{{\mathit{\delta }}_x}}-{\displaystyle \frac{D_{\mathrm{l}}(C_{\mathrm{C}}{\mathit{\mu }}_{\mathrm{C}}-C_{\mathrm{L}}{\mathit{\mu }}_{\mathrm{L}})}{{\mathit{\delta }}_x}}\right).\end{array}\end{array}$$

Again a linearizing assumption $\mathrm{\Delta }\left(\mathit{\mu }C\right)=\partial \left(\mathit{\mu }C\right)/\partial x{\mathit{\delta }}_x$ is made, after which the partial differential equation,

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{\partial \left(C\mathit{\mu }\right)}{\partial t}}={\displaystyle \frac{\partial }{\partial x}}\left(D{\displaystyle \frac{\partial \left(C\mathit{\mu }\right)}{\partial x}}\right),\end{array}$$

is obtained. Deleersnijder et al. (2001) derived this equation with a generalized macroscopic approach.

2.3 Derivation of partial differential equations for higher centralized moments

Centralized moments are defined as

$$\begin{array}{}\text{(8)}& {\mathit{\varphi }}_q={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n({\mathit{\chi }}_i-\mathit{\mu }{)}^q}{n}}.\end{array}$$

The zeroth and first centralized moments are always one and zero, respectively. The variance (σ²), skewness, and other higher moments correspond to q=2, q=3, and q>3. Throughout the text, we shall refer to raw moments as ${\mathit{\mu }}_q=n^{-\mathrm{1}}\sum {\mathit{\chi }}_i^q$ and to non-central moments in general as $n^{-\mathrm{1}}{\sum }_{i=\mathrm{1}}^n({\mathit{\chi }}_i-\mathit{\psi }{)}^q$, where ψ≠μ.

Considering the exchange of matter with the surroundings through the fluxes j_k (Sect. 2.2), the change in q-powered differences in the control volume can be described by

$$\begin{array}{}\text{(9)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{1}}{V}}& \sum _{i=\mathrm{1}}^{n_{\mathrm{n}}}{\left({\mathit{\chi }}_i-{\mathit{\mu }}_{\mathrm{o}}\right)}^q={\displaystyle \frac{\mathrm{1}}{V}}\sum _{j=\mathrm{1}}^{n_{\mathrm{o}}}{\left({\mathit{\chi }}_j-{\mathit{\mu }}_{\mathrm{o}}\right)}^q\\ & +{\displaystyle \frac{\mathrm{1}}{V}}({\mathit{\varphi }}_{j_{\mathrm{1}}}{j}_{\mathrm{1}}-{\mathit{\varphi }}_{j_{\mathrm{2}}}{j}_{\mathrm{2}}-{\mathit{\varphi }}_{j_{\mathrm{3}}}{j}_{\mathrm{3}}+{\mathit{\varphi }}_{j_{\mathrm{4}}}{j}_{\mathrm{4}})\mathrm{\Delta }t,\end{array}\end{array}$$

whereby n_n and n_o denote the number of particles in the updated and old population, respectively. All differences from the mean in Eq. (9), including those associated with mass fluxes, are relative to μ_o. A Taylor series expansion of ϕ around μ_o is used to relate the new state of a population (the left-hand side of Eq. 9) to the new mean age:

$$\begin{array}{}\text{(10)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{1}}{V}}& \sum _{i=\mathrm{1}}^{n_{\mathrm{n}}}{\left({\mathit{\chi }}_i-{\mathit{\mu }}_{\mathrm{n}}\right)}^q={\displaystyle \frac{\mathrm{1}}{V}}\sum _{i=\mathrm{1}}^{n_{\mathrm{n}}}{\left({\mathit{\chi }}_i-{\mathit{\mu }}_{\mathrm{o}}\right)}^q\\ & +C_{\mathrm{n}}{\mathit{\varphi }}^{\prime }\mathrm{\Delta }\mathit{\mu }+C_{\mathrm{n}}{\displaystyle \frac{{\mathit{\varphi }}^{\prime \prime }}{\mathrm{2}}}\mathrm{\Delta }{\mathit{\mu }}^{\mathrm{2}}+C_{\mathrm{n}}{\displaystyle \frac{{\mathit{\varphi }}^{\prime \prime \prime }}{\mathrm{6}}}\mathrm{\Delta }{\mathit{\mu }}^{\mathrm{3}}+\mathrm{\dots },\end{array}\end{array}$$

where $C_{\mathrm{n}}=n_{\mathrm{n}}/V$, $\mathrm{\Delta }\mathit{\mu }={\mathit{\mu }}_{\mathrm{n}}-{\mathit{\mu }}_{\mathrm{o}}$, and ${\mathit{\varphi }}^{\prime }=\partial \mathit{\varphi }/\partial \mathit{\mu }$, etc. The term on the left-hand side of Eq. (10) and the first on the right-hand side of Eq. (9) can be replaced by C_nϕ_n and C_oϕ_o, respectively. By inserting Eq. (10) into Eq. (9) and rearranging the terms, the expression

$$\begin{array}{}\text{(11)}& \begin{array}{rl}{C}_{\mathrm{n}}& {\mathit{\varphi }}_{\mathrm{n}}-C_{\mathrm{o}}{\mathit{\varphi }}_{\mathrm{o}}-C_{\mathrm{n}}{\mathit{\varphi }}^{\prime }\mathrm{\Delta }\mathit{\mu }-C_{\mathrm{n}}{\displaystyle \frac{{\mathit{\varphi }}^{\prime \prime }}{\mathrm{2}}}\mathrm{\Delta }{\mathit{\mu }}^{\mathrm{2}}\\ & -C_{\mathrm{n}}{\displaystyle \frac{{\mathit{\varphi }}^{\prime \prime \prime }}{\mathrm{6}}}\mathrm{\Delta }{\mathit{\mu }}^{\mathrm{3}}-\mathrm{\dots }={\displaystyle \frac{\mathrm{1}}{V}}\left(\sum _{k=\mathrm{1}}^{\mathrm{4}}{\mathit{\lambda }}_k{\mathit{\varphi }}_{j_k}{j}_k\right)\mathrm{\Delta }t\end{array}\end{array}$$

is obtained, whereby ${\mathit{\lambda }}_k=\pm \mathrm{1}$ depending on the direction of the flux.

Table 2 Partial differentials equation for diffusion of concentration, mean, and centralized moments. For the definition of the centralized moments (ϕ_q), see Eq. (8).

2.3.1 Derivation of a partial differential equation for variance

The derivatives of variance in the Taylor series are

$$\begin{array}{}\text{(12a)}& {\displaystyle }& {\displaystyle \frac{\partial {\mathit{\sigma }}^{\mathrm{2}}}{\partial \mathit{\mu }}}=-\mathrm{2}{\displaystyle \frac{\sum ({\mathit{\chi }}_i-\mathit{\mu })}{n}}=\mathrm{0},\text{(12b)}& {\displaystyle }& {\displaystyle \frac{{\partial }^{\mathrm{2}}{\mathit{\sigma }}^{\mathrm{2}}}{\partial {\mathit{\mu }}^{\mathrm{2}}}}=\mathrm{2},\text{(12c)}& {\displaystyle }& {\displaystyle \frac{{\partial }^{\mathrm{3}}{\mathit{\sigma }}^{\mathrm{2}}}{\partial {\mathit{\mu }}^{\mathrm{3}}}}=\mathrm{0}.\end{array}$$

The only non-zero derivative is inserted into Eq. (11). The linearization $\mathrm{\Delta }\mathit{\mu }=\partial \mathit{\mu }/\partial t\mathrm{\Delta }t$ is made, and the result is divided by Δt. Taking the limit of Δt to zero yields

$$\begin{array}{}\text{(13)}& \underset{\mathrm{\Delta }t\to \mathrm{0}}{lim}\left[{\displaystyle \frac{C_{\mathrm{n}}{\mathit{\sigma }}_{\mathrm{n}}^{\mathrm{2}}-C_{\mathrm{o}}{\mathit{\sigma }}_{\mathrm{o}}^{\mathrm{2}}}{\mathrm{\Delta }t}}-C_{\mathrm{n}}{\left({\displaystyle \frac{\partial \mathit{\mu }}{\partial t}}\right)}^{\mathrm{2}}\mathrm{\Delta }t\right]={\displaystyle \frac{\mathrm{1}}{V}}\sum _{k=\mathrm{1}}^{\mathrm{4}}{\mathit{\lambda }}_k{\mathit{\sigma }}_{j_k}^{\mathrm{2}}{j}_k,\end{array}$$

or, under the assumption that $\partial \mathit{\mu }/\partial t$ is finite,

$$\begin{array}{}\text{(14)}& {\displaystyle \frac{\partial \left(C{\mathit{\sigma }}^{\mathrm{2}}\right)}{\partial t}}={\displaystyle \frac{\mathrm{1}}{V}}\sum _{k=\mathrm{1}}^{\mathrm{4}}{\mathit{\lambda }}_k{\mathit{\sigma }}_{j_k}^{\mathrm{2}}{j}_k\end{array}$$

in differential form.

In the next step, the unknown fluxes on the right-hand side are expressed by known local properties, which can only be done for expected mean values of a large number of random particle jumps. It will be assumed for the partial differential equation for variance, as well as for higher-order moments, that (i) the flux is determined by the average jumping frequency and the number of particles from a source location X, i.e., j_k=fn_X; (ii) that q-powered differences reflect the average differences from the location where the particles are jumping, i.e., ${\mathit{\varphi }}_{j_k}={\mathit{\varphi }}_X$; and (iii) that the properties of particles do not effect the jumping probability, i.e., $\langle j_k{\mathit{\varphi }}_{j_k}\rangle =\langle j_k\rangle \langle {\mathit{\varphi }}_{j_k}\rangle$.

With these assumptions, one can write

$$\begin{array}{}\text{(15)}& \begin{array}{rl}& 〈{\displaystyle \frac{\mathrm{1}}{V}}\sum _{k=\mathrm{1}}^{\mathrm{4}}{\mathit{\lambda }}_k{\mathit{\sigma }}_{j_k}^{\mathrm{2}}{j}_k〉={\displaystyle \frac{f_{\mathrm{r}}}{\mathrm{2}}}\left[{\mathit{\sigma }}_{\mathrm{R}}^{\mathrm{2}}\left({\mathit{\mu }}_{\mathrm{C}}\right)C_{\mathrm{R}}-{\mathit{\sigma }}_{\mathrm{C}}^{\mathrm{2}}\left({\mathit{\mu }}_{\mathrm{C}}\right)C_{\mathrm{C}}\right]\\ & -{\displaystyle \frac{f_{\mathrm{l}}}{\mathrm{2}}}\left[{\mathit{\sigma }}_{\mathrm{C}}^{\mathrm{2}}\left({\mathit{\mu }}_{\mathrm{C}}\right)C_{\mathrm{C}}-{\mathit{\sigma }}_{\mathrm{L}}^{\mathrm{2}}\left({\mathit{\mu }}_{\mathrm{C}}\right)C_{\mathrm{L}}\right].\end{array}\end{array}$$

Using the Taylor series for spatial derivatives instead of temporal derivatives, i.e., $\mathrm{\Delta }\mathit{\mu }={\mathit{\mu }}_{\mathrm{R}}-{\mathit{\mu }}_{\mathrm{C}}$ or μ_C−μ_L, according to Eqs. (12a)–(12c) gives the following results:

$$\begin{array}{}\text{(16a)}& {\displaystyle }{\mathit{\sigma }}_{\mathrm{L}}^{\mathrm{2}}\left({\mathit{\mu }}_{\mathrm{C}}\right)={\mathit{\sigma }}_{\mathrm{L}}^{\mathrm{2}}\left({\mathit{\mu }}_{\mathrm{L}}\right)+({\mathit{\mu }}_{\mathrm{C}}-{\mathit{\mu }}_{\mathrm{L}}{)}^{\mathrm{2}},\text{(16b)}& {\displaystyle }{\mathit{\sigma }}_{\mathrm{R}}^{\mathrm{2}}\left({\mathit{\mu }}_{\mathrm{C}}\right)={\mathit{\sigma }}_{\mathrm{R}}^{\mathrm{2}}\left({\mathit{\mu }}_{\mathrm{R}}\right)+({\mathit{\mu }}_{\mathrm{C}}-{\mathit{\mu }}_{\mathrm{R}}{)}^{\mathrm{2}},\end{array}$$

and substituting these into Eq. (15) yields

$$\begin{array}{}\text{(17)}& \begin{array}{rl}& 〈{\displaystyle \frac{\mathrm{1}}{V}}\sum _{k=\mathrm{1}}^{\mathrm{4}}{\mathit{\lambda }}_k{\mathit{\sigma }}_{j_k}^{\mathrm{2}}{j}_k〉=\\ & {\displaystyle \frac{f_{\mathrm{r}}}{\mathrm{2}}}\left[\left({\mathit{\sigma }}_{\mathrm{R}}^{\mathrm{2}}\right({\mathit{\mu }}_{\mathrm{R}})+({\mathit{\mu }}_{\mathrm{C}}-{\mathit{\mu }}_{\mathrm{R}}{)}^{\mathrm{2}})C_{\mathrm{R}}-{\mathit{\sigma }}_{\mathrm{C}}^{\mathrm{2}}({\mathit{\mu }}_{\mathrm{C}})C_{\mathrm{C}}\right]-\\ & {\displaystyle \frac{f_{\mathrm{l}}}{\mathrm{2}}}\left[{\mathit{\sigma }}_{\mathrm{C}}^{\mathrm{2}}\left({\mathit{\mu }}_{\mathrm{C}}\right)C_{\mathrm{C}}-\left({\mathit{\sigma }}_{\mathrm{L}}^{\mathrm{2}}\right({\mathit{\mu }}_{\mathrm{L}})+({\mathit{\mu }}_{\mathrm{C}}-{\mathit{\mu }}_{\mathrm{L}}{)}^{\mathrm{2}})C_{\mathrm{L}}\right].\end{array}\end{array}$$

Ignoring the terms with derivatives obtained from the Taylor series for the moment, a part of the equation can be isolated

$$\begin{array}{}\text{(18)}& \begin{array}{rl}& {〈{\displaystyle \frac{\mathrm{1}}{V}}\sum _{k=\mathrm{1}}^{\mathrm{4}}{\mathit{\lambda }}_k{\mathit{\sigma }}_{j_k}^{\mathrm{2}}{j}_k〉}^{*}={\displaystyle \frac{f_{\mathrm{r}}}{\mathrm{2}}}\left[{\mathit{\sigma }}_{\mathrm{R}}^{\mathrm{2}}\left({\mathit{\mu }}_{\mathrm{R}}\right)C_{\mathrm{R}}-{\mathit{\sigma }}_{\mathrm{C}}^{\mathrm{2}}\left({\mathit{\mu }}_{\mathrm{C}}\right)C_{\mathrm{C}}\right]\\ & -{\displaystyle \frac{f_{\mathrm{l}}}{\mathrm{2}}}\left[{\mathit{\sigma }}_{\mathrm{C}}^{\mathrm{2}}\left({\mathit{\mu }}_{\mathrm{C}}\right)C_{\mathrm{C}}-{\mathit{\sigma }}_{\mathrm{L}}^{\mathrm{2}}\left({\mathit{\mu }}_{\mathrm{L}}\right)C_{\mathrm{L}}\right],\end{array}\end{array}$$

which is similar to Eq. (5). A linearization of $\partial \left({\mathit{\sigma }}^{\mathrm{2}}C\right)/\partial x$ and repeating the procedure leading from Eqs. (5) to (7) here gives

$$\begin{array}{}\text{(19)}& {〈{\displaystyle \frac{\mathrm{1}}{V}}\sum _{k=\mathrm{1}}^{\mathrm{4}}{\mathit{\lambda }}_k{\mathit{\sigma }}_{j_k}^{\mathrm{2}}{j}_k〉}^{*}={\displaystyle \frac{\partial }{\partial x}}\left(D{\displaystyle \frac{\partial \left(C{\mathit{\sigma }}^{\mathrm{2}}\right)}{\partial x}}\right).\end{array}$$

The remaining terms not accounted for yet are

$$\begin{array}{}\text{(20)}& \begin{array}{rl}& {〈{\displaystyle \frac{\mathrm{1}}{V}}\sum _{k=\mathrm{1}}^{\mathrm{4}}{\mathit{\lambda }}_k{\mathit{\sigma }}_{j_k}^{\mathrm{2}}{j}_k〉}^{**}={\displaystyle \frac{f_{\mathrm{r}}}{\mathrm{2}}}({\mathit{\mu }}_{\mathrm{C}}-{\mathit{\mu }}_{\mathrm{R}}{)}^{\mathrm{2}}{C}_{\mathrm{R}}\\ & +{\displaystyle \frac{f_{\mathrm{l}}}{\mathrm{2}}}({\mathit{\mu }}_{\mathrm{C}}-{\mathit{\mu }}_{\mathrm{L}}{)}^{\mathrm{2}}{C}_{\mathrm{L}},\end{array}\end{array}$$

which can also be written as

$$\begin{array}{}\text{(21)}& \begin{array}{rl}& {〈{\displaystyle \frac{\mathrm{1}}{V}}\sum _{k=\mathrm{1}}^{\mathrm{4}}{\mathit{\lambda }}_k{\mathit{\sigma }}_{j_k}^{\mathrm{2}}{j}_k〉}^{**}={\displaystyle \frac{D_{\mathrm{r}}}{{\mathit{\delta }}_x^{\mathrm{2}}}}{\left({\displaystyle \frac{\partial \mathit{\mu }}{\partial x}}{\mathit{\delta }}_x\right)}^{\mathrm{2}}{C}_{\mathrm{R}}\\ & +{\displaystyle \frac{D_{\mathrm{l}}}{{\mathit{\delta }}_x^{\mathrm{2}}}}{\left({\displaystyle \frac{\partial \mathit{\mu }}{\partial x}}{\mathit{\delta }}_x\right)}^{\mathrm{2}}{C}_{\mathrm{L}},\end{array}\end{array}$$

yielding in the limit of Δx→0

$$\begin{array}{}\text{(22)}& {〈{\displaystyle \frac{\mathrm{1}}{V}}\sum _{k=\mathrm{1}}^{\mathrm{4}}{\mathit{\lambda }}_k{\mathit{\sigma }}_{j_k}^{\mathrm{2}}{j}_k〉}^{**}=\mathrm{2}DC{\left({\displaystyle \frac{\partial \mathit{\mu }}{\partial x}}\right)}^{\mathrm{2}}.\end{array}$$

Therefore,

$$\begin{array}{}\text{(23)}& {\displaystyle \frac{\partial \left(C{\mathit{\sigma }}^{\mathrm{2}}\right)}{\partial x}}={\displaystyle \frac{\partial }{\partial x}}\left(D{\displaystyle \frac{\partial \left(C{\mathit{\sigma }}^{\mathrm{2}}\right)}{\partial x}}\right)+\mathrm{2}DC{\left({\displaystyle \frac{\partial \mathit{\mu }}{\partial x}}\right)}^{\mathrm{2}}\end{array}$$

is the final result describing the effect of diffusion on the centralized variance. As demonstrated in Supplement Sect. S2.1, the PDE for diffusion of the raw variance can be derived from Eq. (23) and matches with the result of Delhez and Deleersnijder (2002), which shows that the additional linearizations made in the derivation do not affect the accuracy.

2.3.2 Derivation of partial equations for skewness and all higher-order moments

For a finite $\partial \mathit{\mu }/\partial t$, dividing Eq. (11) by an infinitesimally small time step will drop the higher-order terms in the Taylor series, leaving

$$\begin{array}{}\text{(24)}& {\displaystyle \frac{\partial \left(C\mathit{\varphi }\right)}{\partial t}}-C{\mathit{\varphi }}^{\prime }{\displaystyle \frac{\partial \mathit{\mu }}{\partial t}}={\displaystyle \frac{\mathrm{1}}{V}}\sum _{k=\mathrm{1}}^{\mathrm{4}}{\mathit{\lambda }}_k{\mathit{\varphi }}_{j_k}{j}_k.\end{array}$$

The presence of a non-zero first-order derivative makes the derivation of the PDEs for higher-order moments different from that of the variance. It can be found in Appendix A and is further analytically validated in the Supplement in Sect. S2.2. Table 2 shows an overview of all diffusion PDEs.

3 Derivation of reaction terms for partial differential equations of moments

Here, we will first give a general mathematical approach to derive reaction terms that can be incorporated in the PDEs for centralized moments (Sect. 3.1). The application of this method to particular kinetic expressions relevant to the test applications in this paper will be demonstrated in Sect. 3.2.

3.1 General derivation of differential terms for reactions

Reactions change the concentration and can also influence the shape of distributions characterized by their mean and higher moments. To evaluate the effect of reactions on central moments, we start with an alternative notation for the definition of central moments:

$$\begin{array}{}\text{(25)}& C{\mathit{\varphi }}_q={\displaystyle \frac{\mathrm{1}}{V}}\sum _{p=\mathrm{0}}^q\left[\left(\genfrac{}{}{0ex}{}{q}{p}\right)(-\mathit{\mu }{)}^p\sum _{i=\mathrm{1}}^n{\mathit{\chi }}_i^{q-p}\right],\end{array}$$

which can be obtained by applying the binomial theorem to Eq. (8). Next, we differentiate using the product rule to obtain

$$\begin{array}{}\text{(26)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{d}\left(C{\mathit{\varphi }}_q\right)}{\mathrm{d}t}}=\\ & \sum _{p=\mathrm{0}}^q\left\{(-\mathrm{1}{)}^p\left(\genfrac{}{}{0ex}{}{q}{p}\right)\left[{\mathit{\mu }}^p{\displaystyle \frac{\mathrm{d}\left(C{\mathit{\mu }}_{q-p}\right)}{\mathrm{d}t}}+{\displaystyle \frac{\mathrm{d}{\mathit{\mu }}^p}{\mathrm{d}t}}C{\mathit{\mu }}_{q-p}\right]\right\},\end{array}\end{array}$$

whereby ${\mathit{\mu }}_x=C^{-\mathrm{1}}\int {\mathit{\chi }}^x{C}_{\mathit{\chi }}\mathrm{d}\mathit{\chi }$ can be identified as the xth raw moment, and $C{\mathit{\mu }}_x=\sum {\mathit{\chi }}^x/V$. When no subscript is given, μ=μ₁ denotes the mean. The task at hand is to find expressions for all terms when only the concentration and moments are known from a previous time step during a simulation and a reaction rate is defined, for which we will assume it follows

$$\begin{array}{}\text{(27)}& R=\int r\left(\mathit{\chi }\right)\mathrm{d}\mathit{\chi }\end{array}$$

as a generic rate expression.

Starting with μ_q−p in the second term (Eq. 26), we note that raw moments (μ_x) can be obtained by transforming central moments (ϕ):

$$\begin{array}{}\text{(28)}& {\mathit{\mu }}_x=\sum _{k=\mathrm{0}}^x\left(\genfrac{}{}{0ex}{}{x}{k}\right){\mathit{\varphi }}_k{\mathit{\mu }}^{x-k}.\end{array}$$

The second term is further worked out by substituting

$$\begin{array}{}\text{(29)}& {\displaystyle \frac{\mathrm{d}{\mathit{\mu }}^p}{\mathrm{d}t}}=p{\mathit{\mu }}^{p-\mathrm{1}}{\displaystyle \frac{\mathrm{d}\mathit{\mu }}{\mathrm{d}t}}.\end{array}$$

Solving for $\mathrm{d}\left(C{\mathit{\varphi }}_q\right)/\mathrm{d}t{|}_R$ (Eq. 26) requires expressions for the terms $\mathrm{d}\left(C{\mathit{\mu }}_{q-p}\right)/\mathrm{d}t{|}_R$ and $\mathrm{d}\mathit{\mu }/\mathrm{d}t{|}_R$ (Eq. 29). The first one is obtained by integrating

$$\begin{array}{}\text{(30)}& {\left.{\displaystyle \frac{\mathrm{d}\left(C{\mathit{\mu }}_{q-p}\right)}{\mathrm{d}t}}\right|}_R=\int {\mathit{\chi }}^{q-p}r\left(\mathit{\chi }\right)\mathrm{d}\mathit{\chi }.\end{array}$$

Following the product rule, this derivative also allows

$$\begin{array}{}\text{(31)}& {\left.{\displaystyle \frac{\mathrm{d}\mathit{\mu }}{\mathrm{d}t}}\right|}_R=\left[{\left.{\displaystyle \frac{\mathrm{d}\left(C\mathit{\mu }\right)}{\mathrm{d}t}}\right|}_R-\mathit{\mu }R\right]C^{-\mathrm{1}}\end{array}$$

to be solved.

3.2 Implementation of reaction kinetics and aging

Representing the production of new material with a zeroth-order kinetic term (r=0) results in R=P as the constant of integration when treating the integral in Eq. (27) as indefinite. When χ represents age, the PDEs for raw moments times concentration do not need additional terms to account for production (${\left.\mathrm{d}\left[C{\mathit{\mu }}_q\right]/\mathrm{d}t\right|}_R=\mathrm{0}$ since χ=0, Eq. 30). However, the production will decrease the mean age (Eq. 31). To incorporate production into the PDEs of variance and higher central moments, the additional term P_q is obtained by inserting ${\left.\mathrm{d}\left(C{\mathit{\mu }}_q\right)/\mathrm{d}t\right|}_R=\mathrm{0}$ and $\mathrm{d}\mathit{\mu }/\mathrm{d}t{|}_R=-\mathit{\mu }P{C}^{-\mathrm{1}}$ into Eq. (26).

The formulation

$$\begin{array}{}\text{(32)}& r=f\left(\mathit{\chi }\right)C\left(\mathit{\chi }\right)\end{array}$$

describes first-order reaction kinetics, whereby f(χ) is the reactivity as a function of the distributed property, and the integrations in Eqs. (27) and (30) are performed over the domain bounded by the scope of the χ distribution. However, when $f\left(\mathit{\chi }\right)=-k$ is constant, the reaction rate only depends on concentration ($R=-kC$), as is the case for radioactive decay and simple first-order kinetics. As these reactions do not discriminate with respect to age, the moments will not change (i.e., $\mathrm{d}{\mathit{\varphi }}_q/\mathrm{d}t{|}_R=\mathrm{d}{\mathit{\mu }}_q/\mathrm{d}t{|}_R=\mathrm{0}$). Consequently, the reaction term becomes $\mathrm{d}\left(C{\mathit{\varphi }}_q\right)/\mathrm{d}t{|}_R=R{\mathit{\varphi }}_q$, following the product rule.

The distribution may directly represent the reactivity, e.g., $f\left(\mathit{\chi }\right)=-\mathit{\chi }$ (see the application in Sect. 4.2). When in total x moments are simulated, functions of χ^x+1 will have to be integrated (Eqs. 30, 32). For simulations whereby the integrations are numerically performed, this may substantially impact the computation time.

In the final application (Sect. 4.3), a hypothetical first-order reaction rate will be considered, whereby the reactivity function (f[χ]) depends on the inverse of χ, which will represent age. Aging affects the mean but not the concentration. Hence, the product rule implies $\mathrm{d}\left(C\mathit{\mu }\right)/\mathrm{d}t=C\mathrm{d}\mathit{\mu }/\mathrm{d}t$, whereby $\mathrm{d}\mathit{\mu }/\mathrm{d}t$ due to aging will be unity provided the same units are used for χ and t. Aging shifts the distribution along the age axis but does not change the shape of the distribution. It does not contribute to changes in the differences between particle ages and the mean ($\mathrm{d}({\mathit{\chi }}_i-\mathit{\mu }{)}^q=\mathrm{0}$), nor does it impact the concentration. This implies that $\mathrm{d}\left(C{\mathit{\varphi }}_q\right)/\mathrm{d}t$ due to aging is zero.

4 Applications

Three applications related to sedimentary environments are presented. For the sake of simplicity and generality, the effect of sediment properties, such as porosity and tortuosity, on transport will be ignored. Instead, the focus is on adding reactions.

4.1 Simulating an age tracer

In sediment modeling, resolving the age of a chemical compound understood as the time since its formation or deposition onto the sediment can help to fit a measured profile and serve other diagnostic purposes. A general system of equations for an age associated with a chemical concentration C is given by

$$\begin{array}{}\text{(33a)}& {\displaystyle \frac{\partial C}{\partial t}}& {\displaystyle }=J_{\mathrm{0}}-{\displaystyle \frac{\partial \left(\mathit{\omega }C\right)}{\partial x}}+P+R,\text{(33b)}& {\displaystyle \frac{\partial \left(C\mathit{\mu }\right)}{\partial t}}& {\displaystyle }=J_{\mathrm{1}}-{\displaystyle \frac{\partial \left(\mathit{\omega }C\mathit{\mu }\right)}{\partial x}}+C+\mathit{\mu }R,\text{(33c)}& {\displaystyle \frac{\partial \left(C{\mathit{\varphi }}_q\right)}{\partial t}}& {\displaystyle }=J_q-{\displaystyle \frac{\partial \left(\mathit{\omega }C{\mathit{\varphi }}_q\right)}{\partial x}}+P_q+{\mathit{\varphi }}_{q}R.\end{array}$$

In the equations, J denotes the diffusive transport terms listed in Table 2. The second term accounts for advective transport, whereby ω is the velocity. In early diagenetic models, the accumulation of the sediment column is typically described as a downward advective burial process, since the sediment surface stays at a zero vertical coordinate. The term P denotes the production of new material with an age of zero, P_q accounts for the effect of production on higher centralized moments, R represents a consumption reaction that does not discriminate with respect to age (e.g., $R=-kC$), and the third term in Eq. (33b) accounts for aging. Refer to Sect. 3.2 for the derivation of these terms.

Figure 1 Profiles of moments evolving over time due to diffusion. Solution of numerical integration of partial differential equations (solid lines) is compared to a particle tracking simulation (circles). Domain contains 50 cells. The spacing and time steps are set equal to the jumping distance and the inverse jumping frequency of particles, respectively.

An example of a simulation involving diffusion without aging, advection, and reactions is shown in Fig. 1. Here fixed concentrations and moments were imposed for the last and first cells in the domain as boundary conditions. The initial conditions for the domain were set to the left boundary condition, which has a lower concentration and different moments compared to the right boundary. Over time, there is net chemical diffusion in the leftward direction throughout the domain, eventually leading to a new steady state. The computed concentration and the first seven moments match well with those computed by a particle-based simulation. There is a small but noticeable mismatch at the peak for the skewness and higher moments, which is potentially due to the finite step size in the Lagrangian simulation.

The Eulerian simulation is based on a finite-differences scheme, implemented in R (R Core Team, 2022) and run with the CVODE solver (Brown et al., 1989; Soetaert et al., 2010). The Lagrangian model employed for validating the Eulerian simulation is described in the Supplement Sect. S1.1. The script to run these simulations is relatively simple and publicly available online. Therein, it is possible to add reactions for production and consumption.

4.2 Simulating organic matter mineralization with a reactivity continuum model in turbated sediments

In this application, C denotes organic carbon concentration, χ is the reactivity (degradation rate coefficient) with dimensions T⁻¹, and no explicit aging process is involved. As an example, we will consider the deposition of organic carbon with an initially uniform distribution for $\mathit{\chi }\in [\mathrm{0},m]$, which is described by the state variables concentration, mean reactivity, and variance of the reactivity. The rate expression from Eq. (32) is applied with $f\left(\mathit{\chi }\right)=-\mathit{\chi }$ so that $R=-\mathit{\mu }C$. By working out Eq. (26) with these definitions, the following equations can be obtained.

$$\begin{array}{}\text{(34a)}& {\displaystyle }{\left.{\displaystyle \frac{\partial C}{\partial t}}\right|}_R=-\mathit{\mu }C\text{(34b)}& {\displaystyle }{\left.{\displaystyle \frac{\partial \left(C\mathit{\mu }\right)}{\partial t}}\right|}_R=-\underset{\mathrm{0}}{\overset{m}{\int }}\left({\mathit{\chi }}^{\mathrm{2}}C\right)\mathrm{d}\mathit{\chi }\text{(34c)}& {\displaystyle }{\left.{\displaystyle \frac{\partial \left(C{\mathit{\sigma }}^{\mathrm{2}}\right)}{\partial t}}\right|}_R=-\underset{\mathrm{0}}{\overset{m}{\int }}\left({\mathit{\chi }}^{\mathrm{3}}C\right)\mathrm{d}\mathit{\chi }+\mathrm{2}\mathit{\mu }\underset{\mathrm{0}}{\overset{m}{\int }}\left({\mathit{\chi }}^{\mathrm{2}}C\right)\mathrm{d}\mathit{\chi }-{\mathit{\mu }}^{\mathrm{3}}C\end{array}$$

The integrals will be evaluated numerically, meaning that the full distribution needs to be constructed from the moments. Based on a finite number of moments, it can only be estimated, and we chose the following function to represent the reactivity distribution.

$$\begin{array}{}\text{(35)}& g(\mathit{\chi },\mathit{w})=C_{\mathrm{0}}{e}^{w_{\mathrm{0}}\sqrt{\mathit{\chi }}+w_{\mathrm{1}}\mathit{\chi }}+w_{\mathrm{2}}\mathit{\chi }{e}^{\mathit{\chi }}\end{array}$$

It is motivated as follows: the concentration of unreactive organic matter at the intercept does not change ($g(\mathrm{0},\mathit{w})=C_{\mathrm{0}}$). For the diffusion–reaction equation $\partial C/\partial t=D{\partial }^{\mathrm{2}}C/\partial x^{\mathrm{2}}-\mathit{\chi }C$, the general solution is $A{e}^{\pm \sqrt{\mathit{\chi }/D}x}$. The solution for $\partial C/\partial t=-\mathit{\chi }C$ is $C\left(t\right)=A{e}^{-\mathit{\chi }t}$. The first term can capture both these dynamics. The last term is linearly independent and introduces a third fitting parameter to match the number of equations. These terms have the desirable properties that they cannot fluctuate or become negative and can be evaluated at χ=0. The following equations are solved with a multidimensional root-finding procedure (Soetaert, 2009) to fit the parameter vector w.

$$\begin{array}{}\text{(36a)}& {\displaystyle }C-\underset{\mathrm{0}}{\overset{m}{\int }}g(\mathit{\chi },\mathit{w})\mathrm{d}\mathit{\chi }=\mathrm{0}\text{(36b)}& {\displaystyle }\mathit{\mu }-\underset{\mathrm{0}}{\overset{m}{\int }}g(\mathit{\chi },\mathit{w})\mathit{\chi }\mathrm{d}\mathit{\chi }=\mathrm{0}\text{(36c)}& {\displaystyle }{\mathit{\sigma }}^{\mathrm{2}}-\underset{\mathrm{0}}{\overset{m}{\int }}(\mathit{\chi }-\mathit{\mu }{)}^{\mathrm{2}}g(\mathit{\chi },\mathit{w})\mathrm{d}\mathit{\chi }=\mathrm{0}\end{array}$$

Figure 2 A simulated reactivity continuum model run to a steady state (black line), validated against a model using 30 bins to encompass the reactivity range (dashed light blue line). Concentration, mean reactivity, and standard deviation are shown in panels (a), (b), and (c), respectively. Panel (d) depicts the distribution function for reactivities at the upper (green line) and lower (brown line) boundaries of the model domain in comparison to the validation simulation (points).

In the example shown in Fig. 2, transport involves bioturbation and advection. The burial velocity was set to 1 mm yr⁻¹. The bioturbation coefficient had a maximum value of 10⁻¹⁰ m² s⁻¹ at the sediment–water interface and decreased exponentially over depth with an e-folding distance set to 2 cm. The uniform distribution implemented as Dirichlet upper-boundary conditions was defined by the moments C=1, $\mathit{\mu }=\mathrm{2.5}\times {\mathrm{10}}^{-\mathrm{2}}$ yr⁻¹, and $\mathit{\sigma }=\mathrm{1.44}\times {\mathrm{10}}^{-\mathrm{2}}$ yr⁻¹, which were also used as initial conditions throughout the domain. A no-gradient condition was set as lower-boundary condition. In total, 50 evenly spaced cells discretized a domain length of 50 cm. The simulation, using finite differences (Soetaert and Meysman, 2012), was run with the VODE solver (Brown et al., 1989).

The organic matter concentration imposed to a fixed value at the upper-boundary condition decreases due to degradation over depth (Fig. 2a). Due to mixing, the mean reactivity of organic matter remains relatively constant in the bioturbated zone but decreases below it (Fig. 2b) as the more reactive organic matter is degraded. The variance is also kept relatively stable within the bioturbated zone and decreases strongly below (Fig. 2c), as the removal of more reactive material decreases the spreading of the reactivity distribution. The distributions at the top and bottom are shown as well (Fig. 2d). The obtained results closely resemble those of a 30-G model, which partitions the reactivity range defined at the upper boundary into 30 equally spaced distinct reactivity values, treating each bin as an independent state variable (see for more details Sect. S1.2 of the Supplement).

4.3 Apparent organic matter reactivity as a function of age

In this application, the age distribution is modeled and determines the reactivity of organic carbon. The transport equations are the same as in the previous applications, and aging is included (see Eq. 33b). The age-dependent reactivity (Eq. 32) is specified as

$$\begin{array}{}\text{(37)}& f\left(\mathit{\chi }\right)=-{\displaystyle \frac{v}{\mathit{\alpha }+\mathit{\chi }}},\end{array}$$

whereby α prevents division by zero and may also be used as a fitting parameter. It resembles the long-established expression for the mean reactivity, $\overline{k}=v/(\mathit{\alpha }+\mathit{\mu }{)}^{\mathit{\beta }}$ (Middelburg, 1989). Conceptually it may be the most simple expression to model the effect of aging on reactivity.

Here we consider the application of four moments: concentration (C), mean age (μ), variance (σ²=ϕ₂), and skewness (S=ϕ₃). For this, three-parameter distributions describe the distribution's shape and mean, and a fourth parameter serves as a multiplier to adjust the concentration.

We present our analysis using two distinct distributions to represent age continua, which will be compared later to evaluate the distribution shape's role. The first one is the triangular distribution:

$$\begin{array}{}\text{(38)}& g(\mathit{\chi },\mathit{w})=\left\{\begin{array}{ll}{w}_{\mathrm{1}}\frac{\mathrm{2}(\mathit{\chi }-w_{\mathrm{2}})}{(w_{\mathrm{3}}-w_{\mathrm{2}})(w_{\mathrm{4}}-w_{\mathrm{2}})}& \text{if }{w}_{\mathrm{2}}\le \mathit{\chi }\le w_{\mathrm{4}}\\ w_{\mathrm{1}}\frac{\mathrm{2}(w_{\mathrm{3}}-\mathit{\chi })}{(w_{\mathrm{3}}-w_{\mathrm{2}})(w_{\mathrm{3}}-w_{\mathrm{4}})}& \text{if }{w}_{\mathrm{4}}<\mathit{\chi }\le w_{\mathrm{3}}\end{array}\right.,\end{array}$$

whereby w₁ is the multiplier, w₂ and w₃ denote the lower and upper limit of the distribution (outside this interval, the function evaluates to zero), and w₄ corresponds to the mode. Given the closed-form expressions for the central moments of the triangular distribution (Forbes et al., 2010), the distribution parameters for given moments are found by first determining b and c in the following equations with a root solver

$$\begin{array}{}\text{(39a)}& {\displaystyle }{\mathit{\sigma }}^{\mathrm{2}}-(b^{\mathrm{2}}+c^{\mathrm{2}}-bc)/\mathrm{18}=\mathrm{0},\text{(39b)}& {\displaystyle }S+(b-\mathrm{2}c)(c-\mathrm{2}b)(b+c)/\mathrm{270}=\mathrm{0},\end{array}$$

which allows the distribution parameters to be calculated as follows: w₁=C, $w_{\mathrm{2}}=(\mathrm{3}\mathit{\mu }-b-c)/\mathrm{3}$, $w_{\mathrm{3}}=w_{\mathrm{2}}+b$, and $w_{\mathrm{4}}=w_{\mathrm{2}}+c$.

The other demonstrated distribution is the translated Weibull distribution, formulated as

$$\begin{array}{}\text{(40)}& g(\mathit{\chi },\mathit{w})={\displaystyle \frac{w_{\mathrm{1}}{w}_{\mathrm{2}}}{w_{\mathrm{3}}}}{\left({\displaystyle \frac{\mathit{\chi }-w_{\mathrm{4}}}{w_{\mathrm{3}}}}\right)}^{w_{\mathrm{2}}-\mathrm{1}}{e}^{-{\left(\frac{\mathit{\chi }-w_{\mathrm{4}}}{w_{\mathrm{3}}}\right)}^{w_{\mathrm{2}}}},\end{array}$$

whereby w₁ serves as a multiplier to adjust the concentration, w₂ is the shape parameter, w₃ is the scaling parameter, and w₄ is the location parameter (Forbes et al., 2010). All parameters were obtained by solving a set of equations as shown in Eq. (36), along with an additional equation to account for skewness.

The numerical model first calculates the reaction rates for all moments, as described in Sect. 3.2. The PDEs are solved with an implicit finite-volume scheme with hybrid differences to account for advection and diffusion using the implementation from JurRTM (Rooze et al., 2020; Zindorf et al., 2021). For the state variable Cμ, the aging term is added to the reaction term. In addition, the last term in the diffusion equations for variance and skewness (Table 2) is accounted for by calculating first the $\partial \mathit{\mu }/\partial x$ gradient and adding the resulting term as a reaction rate.

The model divides a domain length of 10 cm into 50 evenly spaced cells. The upper-boundary condition is added as a Dirichlet boundary condition with a prescribed distribution. The initial values are set to the values of the upper-boundary condition. A zero-gradient condition is imposed at the lower boundary of the domain. The reactivity parameters (Eq. 37) were set to v=0.8 and α=1 year.

The simulations were run for 53 years with a maximum time step of 1 year. For root solving the “nleqslv” and “numDeriv” packages in R were utilized (Hasselman, 2023; Gilbert and Varadhan, 2016; R Core Team, 2022). Numerical integrations were carried out over the domain $\mathit{\chi }\in [\mathrm{0},\mathrm{150}]$. To validate the simulation, a complementary model capturing concentration evolution across binned ages was employed. This validation model adopts matching boundary and initial conditions, along with a congruent simulation setup. Detailed technical documentation for the validation model is provided in the Supplement (Sect. S1.3).

Figure 3 Results of simulations based on age distributions. Panels (a) and (b) depict simulations utilizing the triangular and translated Weibull distributions, respectively. Presented are results after both 15 and 50 years of simulation. The moments-based simulation is compared to a simulation using age bins. The moment-based and age bin simulations are distinguished by solid and dashed lines in the upper two rows (see legend) and by lines and points in the lower two rows, which display distribution functions (df). For additional details on the simulations, refer to Sect. 4.3.

For the simulations shown in Fig. 3, the burial velocity was set to 2 mm yr⁻¹. The bioturbation coefficient had a maximum value of 10⁻¹¹ m² s⁻¹ at the top. It decreased exponentially with depth, having an e-folding distance set to 3 cm, implying that the diffusivity at the bottom of the domain is effectively zero and that the zero-gradient condition has a negligible effect on the results. The distributions set as upper-boundary and initial conditions had moments set to C=1, μ=15 years, σ²=13 yr², and S=0 yr³.

The results (Fig. 3) demonstrate accurate concentration and mean age computation throughout the simulations. However, a noticeable mismatch appears in the simulated variance and skewness after 50 years in the triangular distribution simulation, whereas the translated Weibull distribution simulation reproduces the moments accurately.

In the two lower rows of Fig. 3, distributions are presented as resolved during simulations for integration. After 15 years, the moment-based simulation distributions closely match the validation simulation, except at 2.5 cm depth for the triangular distribution. Mixing of younger and older materials from the upper boundary results in a positive skew near the upper boundary and a negative skew at greater depths.

The influence of the initial conditions diminishes over time. Considering only advection, the time to transport all initial material out of the model domain is 50 years. However, upward diffusion, despite net downward chemical diffusion, can increase the residence time of some particles. The triangular distribution struggles to accurately depict the resulting positively skewed distributions, as it degenerates into nearly right triangles at and below 2.5 cm depth (Fig. 3). The much more versatile Weibull distribution provides a better representation of the age distribution, but the distribution is clearly off in the most actively turbated zone (Fig. 3), which, however, does not appear to affect the accuracy of the simulated moments. When the simulation is run for 100 years (not shown), the moments maintain their accuracy, and the visual comparison of the distribution even slightly improves.

Interestingly, the limited impact of the precise shape of the distribution also transpires from the great similarity of the distributions that evolve at depth over time in the validation simulations, regardless of the imposed distribution type at the upper-boundary condition (compare distributions below 2.5 cm depth after 50 years emerging from triangular and Weibull distributions in Fig. 3). Similar distributions also formed when the reaction was turned off (not shown). Therefore, the interplay between bioturbation and aging appeared to be of greater consequence for the evolution of the moments than the reaction and choice of distribution at the upper boundary.

5 Discussion

5.1 Evaluation of the applications

The theory outlined in this paper gives modelers a free hand to simulate properties associated with particles described as a concentration. Modelers can use different reaction kinetics and simulate the effect on a continuous distribution. The first application shows that it is possible to reliably simulate numerous moments of an age distribution. In this type of application, it is unnecessary to reconstruct the distribution function from its moments during the simulation. In principle, any distribution is uniquely defined by a large or infinite number of moments. However, in practice, there is no universal solution to retrieve distributions from their moments. Numerical methods have been developed for this purpose (e.g., John et al., 2007; Arbel et al., 2016) but do not always succeed.

In the two other applications, the shape of the distribution affected the reaction rates, and it therefore needed to be resolved in the simulations. In the reactivity continuum model, the shape could be well predicted; i.e., the chosen class of distribution function was a good approximation of the exact solutions. The bounds of the initial distribution of deposited material also provided bounds for the numerical integration of deeper (older) material since the domain only becomes smaller as the more reactive materials are consumed, and the refractory organics remain. These simulations were stable and ran relatively fast, despite the requirement to determine the parameters of the distribution in runtime, which involves several numerical integrations. This approach could be attractive as an alternative to the multi-G approach, as it does not require the arbitrary definition of various reactivity classes (Jørgensen, 1978) and can better represent slight differences in reactivity naturally developing over depth. However, the multi-G simulation will run faster, even though the numerical scheme adopted here could be significantly improved (see suggestions below). The initial uniform distribution used in the simulation may not be realistic for organic matter in marine sediment (Boudreau and Ruddick, 1991). Instead, distributions similar to multi-G or continuous distributions found in the literature (Arndt et al., 2013) can be imposed for freshly deposited organics.

The third application posed the greatest challenge, as the shape of the age distribution strongly changed during the simulation and was hard to predict. The emergence of pronounced asymmetry motivated the addition of skewness as a state variable. The translated Weibull distribution performed best due to its versatility, as it can represent the exponential, Rayleigh, normal, and other two-parameter distributions.

Occasional substantial deviations were observed between distributions reconstructed from their moments and the actual distributions. An illustration of such discrepancies is evident in the left tails of the distributions after 50 years at a depth of 2.5 cm (Fig. 3), even though the moments of both the triangular and translated Weibull distributions were nearly or entirely accurate. A fundamental problem of reconstructing distributions from moments pertains to the unequal weighting for concentration differences along the age–reactivity axis. Central moments exhibit zero weight at χ=μ (Eq. 8), while weights are maximized at the distribution's extremes. Consequently, the reconstruction process tends to fit the tails better than the central region encompassing the mean and mode. This tendency could introduce biases, particularly in scenarios where the extremes influence the dynamics less; for instance, when the precise age of very old material has minimal impact on overall reactivity. The bias, depending on q in Eq. (8), will be stronger for higher-order moments.

Numerous optimizations can be considered to improve the numerical schemes. In the third simulation, a transition from the method of lines (Sarmin and Chudov, 1963), characterized by continuous time and discrete space, to the classical finite-differences and volumes approach, characterized by discrete space and time, yielded substantial improvements in simulation time and stability. This approach, affording greater control over time stepping and the execution frequency of root-solving procedures and numerical integrations could likely also shorten computation times in the other applications. Other polynomial type or spline expressions could be tried to describe the distributions. Finally, pre-calculated search tables for distribution functions could be designed to look up parameters corresponding to a combination of moments. Then numerical integrations during run time would become obsolete, letting simulations run faster.

5.2 The application of central moments-based models in comparison to alternative approaches

Instead of utilizing central moments, an alternative consideration involves using raw moments. When focusing solely on production processes and age distributions are simulated, raw moments could prove more practical, as they are not affected by the production of new material, but the disadvantage will be that aging will still affect them.

In scenarios involving consumption reactions, the use of raw moments generally leads to the evaluation of the fewest number of terms. For example, consider the rate expression in Eq. (30) compared to those in Eq. (34) for central moments. However, the steps outlined in Sect. 3.1 can be automated to obtain all necessary terms. For the complete set of equations encompassing all moments, the same integrations must be carried out, regardless of whether raw or central moments are employed. Hence, the choice between raw and central moments may have limited practical significance when considering only the consumption reactions in numerical models.

Central moments can be converted into raw moments (as in Eq. 28) or any other non-central moment. This implies that the choice of moment type for the PDEs is not critical for additional steps, such as the reconstruction of distributions from the moments (see previous section).

Central moments have advantages in simulating age or transit time distributions, as it lets aging only affect the mean and not higher moments like variance and skewness. In contrast, PDEs for non-central moments necessitate additional terms to account for aging (Delhez and Deleersnijder, 2002). While moments-based distributions have been employed for simulating age tracers and radioactive decay, their application in more complex dynamics remains unexplored. Determining the practical benefits of central versus non-central moments would require further testing.

For each application within this study, alternative numerical approaches were presented. One such approach involves utilizing Lagrangian simulation, which typically demands more computation time and may be less suitable for boundary value problems involving extended simulation durations needed to reach a steady state. In the second and third applications, continuous distributions were discretized by multiple state variables. While multi-G models run faster and are easier to implement, continuous approaches become particularly relevant when the goal is to examine aging processes. Discretizing age distributions presents challenges, as aging involves exchanges between different age bins. Similar to numerical advection schemes, this can lead to numerical diffusion, distorting the variance and skewness (Klingbeil et al., 2014). To circumvent this concern, the validation simulation in the third application utilized a moving grid for age bins (Supplement Sect. S1.3). This approach reaches high accuracy. However, it has the disadvantage of requiring many state variables to represent age classes for the simulated period, which could become problematic, particularly in simulations with larger grids. Moments-based simulations offer an elegant and efficient solution, while the versatile applicability of PDEs makes their implementation in various models more convenient.

6 Conclusions

The derived diffusion–advection–reaction PDEs for central moments can be valuable tools for assessing the effect of processes on distributions, computing transit time and age distributions, and simulating reactivity distributions. Central moments hold advantages over raw moments, being intuitively interpretable and unaffected by aging.

The central moments of transit time and age distributions can be simulated first to allow the actual distribution reconstruction afterward. When the reactivity depends on the distribution, the distribution must be reconstructed and integrated at each time step. An adequate function could be defined to carry out the reconstruction from the mean and variance for the simulated distributions representing the reactivity continua in the second application. However, resolving age distributions in the third application to compute reactivities based on ages encountered distribution-choice-related accuracy challenges, suggesting a need for additional validation, particularly when the distribution is more sensitive to the reaction term and vice versa.

The second and third applications underscored bioturbation's considerable influence on chemical reactivity and age distributions in surficial sediments, highlighting potential inaccuracies in prior reactivity continuum approaches that ignore mixing. Despite employing realistic transport parameters, applying the models to field data remains essential, particularly for more robust validation of the chosen reaction dynamics. In addition, a more thorough analysis is required to assess the significance of age and reactivity distribution shapes for mineralization rates within mixed zones. The framework developed within this study is well suited to address these aspects in future research.

Appendix A: Derivation of diffusion PDEs for higher centralized moments

Continuing the derivation of the PDE for the diffusion of higher moments from Eq. (24), one can write the following equation:

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{1}}{V}}〈\sum {\mathit{\lambda }}_k{\mathit{\varphi }}_{j_k}{j}_k〉& =-{\displaystyle \frac{f_{\mathrm{l}}}{\mathrm{2}}}\left[{\mathit{\varphi }}_{\mathrm{C}}\left({\mathit{\mu }}_{\mathrm{C}}\right)C_{\mathrm{C}}-{\mathit{\varphi }}_{\mathrm{L}}\left({\mathit{\mu }}_{\mathrm{C}}\right)C_{\mathrm{L}}\right]\\ & +{\displaystyle \frac{f_{\mathrm{r}}}{\mathrm{2}}}\left[{\mathit{\varphi }}_{\mathrm{R}}\left({\mathit{\mu }}_{\mathrm{C}}\right)C_{\mathrm{R}}-{\mathit{\varphi }}_{\mathrm{C}}\left({\mathit{\mu }}_{\mathrm{C}}\right)C_{\mathrm{C}}\right],\end{array}\end{array}$$

where the same assumptions are made with regard to the fluxes as in the derivation for the PDE of the variance. The fluxes j₁ and j₄, transporting material into the control volume, can be written as functions of the mean age at the source location,

$$\begin{array}{}\text{(A2a)}& {\displaystyle }\begin{array}{rl}{\mathit{\varphi }}_{\mathrm{L}}\left({\mathit{\mu }}_{\mathrm{L}}\right)& ={\mathit{\varphi }}_{\mathrm{L}}\left({\mathit{\mu }}_{\mathrm{C}}\right)+{\mathit{\varphi }}^{\prime }\cdot ({\mathit{\mu }}_{\mathrm{L}}-{\mathit{\mu }}_{\mathrm{C}})+\mathrm{0.5}{\mathit{\varphi }}^{\prime \prime }\\ & \cdot ({\mathit{\mu }}_{\mathrm{L}}-{\mathit{\mu }}_{\mathrm{C}}{)}^{\mathrm{2}}+\mathrm{\dots },\end{array}\text{(A2b)}& {\displaystyle }\begin{array}{rl}{\mathit{\varphi }}_{\mathrm{R}}\left({\mathit{\mu }}_{\mathrm{R}}\right)& ={\mathit{\varphi }}_{\mathrm{R}}\left({\mathit{\mu }}_{\mathrm{C}}\right)+{\mathit{\varphi }}^{\prime }\cdot ({\mathit{\mu }}_{\mathrm{R}}-{\mathit{\mu }}_{\mathrm{C}})+\mathrm{0.5}{\mathit{\varphi }}^{\prime \prime }\\ & \cdot ({\mathit{\mu }}_{\mathrm{R}}-{\mathit{\mu }}_{\mathrm{C}}{)}^{\mathrm{2}}+\mathrm{\dots }.\end{array}\end{array}$$

Inserting these equations into Eq. (A1), the part not accounting for the derivatives of the Taylor series is isolated as follows:

$$\begin{array}{}\text{(A3)}& {\displaystyle \frac{\mathrm{1}}{V}}{〈\sum {\mathit{\lambda }}_k{\mathit{\varphi }}_{j_k}{j}_k〉}^{*}={\displaystyle \frac{\partial }{\partial x}}\left(D{\displaystyle \frac{\partial \left(C\mathit{\varphi }\right)}{\partial x}}\right),\end{array}$$

whereby $\partial \left(C\mathit{\varphi }\right)/\partial x$ has been linearized. The terms for the derivatives can be written as follows:

$$\begin{array}{}\text{(A4)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{1}}{V}}{〈\sum {\mathit{\lambda }}_k{\mathit{\varphi }}_{j_k}{j}_k〉}^{**}=\\ & -{\displaystyle \frac{f_{\mathrm{l}}{C}_{\mathrm{L}}}{\mathrm{2}}}\left({{\displaystyle \frac{\partial \mathit{\varphi }}{\partial \mathit{\mu }}}}_{\mathrm{L}}({\mathit{\mu }}_{\mathrm{L}}-{\mathit{\mu }}_{\mathrm{C}})+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\varphi }}{\partial {\mathit{\mu }}^{\mathrm{2}}}}}_{\mathrm{L}}({\mathit{\mu }}_{\mathrm{L}}-{\mathit{\mu }}_{\mathrm{C}}{)}^{\mathrm{2}}+\mathrm{\dots }\right)-\\ & {\displaystyle \frac{f_{\mathrm{r}}{C}_{\mathrm{R}}}{\mathrm{2}}}\left({{\displaystyle \frac{\partial \mathit{\varphi }}{\partial \mathit{\mu }}}}_{\mathrm{R}}({\mathit{\mu }}_{\mathrm{R}}-{\mathit{\mu }}_{\mathrm{C}})+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\varphi }}{\partial {\mathit{\mu }}^{\mathrm{2}}}}}_{\mathrm{R}}({\mathit{\mu }}_{\mathrm{R}}-{\mathit{\mu }}_{\mathrm{C}}{)}^{\mathrm{2}}+\mathrm{\dots }\right).\end{array}\end{array}$$

Substituting this and $f/\mathrm{2}=D/{\mathit{\delta }}_x^{\mathrm{2}}$ into the last equation yields

$$\begin{array}{}\text{(A5)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{1}}{V}}{〈\sum {\mathit{\lambda }}_k{\mathit{\varphi }}_{j_k}{j}_k〉}^{**}={\displaystyle \frac{D_{\mathrm{l}}{C}_{\mathrm{L}}}{{\mathit{\delta }}_x}}\\ & \left({{\displaystyle \frac{\partial \mathit{\varphi }}{\partial \mathit{\mu }}}}_{\mathrm{L}}{{\displaystyle \frac{\mathrm{\Delta }\mathit{\mu }}{{\mathit{\delta }}_x}}}_{\mathrm{l}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\varphi }}{\partial {\mathit{\mu }}^{\mathrm{2}}}}}_{\mathrm{L}}{{\displaystyle \frac{(\mathrm{\Delta }\mathit{\mu }{)}^{\mathrm{2}}}{{\mathit{\delta }}_x}}}_{\mathrm{l}}+\mathrm{\dots }\right)\\ & -{\displaystyle \frac{D_{\mathrm{r}}{C}_{\mathrm{R}}}{{\mathit{\delta }}_x}}\left({{\displaystyle \frac{\partial \mathit{\varphi }}{\partial \mathit{\mu }}}}_{\mathrm{R}}{{\displaystyle \frac{\mathrm{\Delta }\mathit{\mu }}{{\mathit{\delta }}_x}}}_{\mathrm{r}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\varphi }}{\partial {\mathit{\mu }}^{\mathrm{2}}}}}_{\mathrm{R}}{{\displaystyle \frac{(\mathrm{\Delta }\mathit{\mu }{)}^{\mathrm{2}}}{{\mathit{\delta }}_x}}}_{\mathrm{r}}+\mathrm{\dots }\right),\end{array}\end{array}$$

whereby ${\mathit{\mu }}_{\mathrm{C}}-{\mathit{\mu }}_{\mathrm{L}}=\mathrm{\Delta }{\mathit{\mu }}_{\mathrm{l}}$ and ${\mathit{\mu }}_{\mathrm{R}}-{\mathit{\mu }}_{\mathrm{C}}=\mathrm{\Delta }{\mathit{\mu }}_{\mathrm{r}}$. Taking the limit of Δx to zero, the second-order Taylor series terms will drop. Linearizing $\partial \mathit{\mu }/\partial x$ yields

$$\begin{array}{}\text{(A6)}& {\displaystyle \frac{\mathrm{1}}{V}}{〈\sum {\mathit{\lambda }}_k{\mathit{\varphi }}_{j_k}{j}_k〉}^{**}={\displaystyle \frac{D_{\mathrm{l}}{C}_{\mathrm{L}}}{{\mathit{\delta }}_x}}{{\displaystyle \frac{\partial \mathit{\varphi }}{\partial \mathit{\mu }}}}_{\mathrm{L}}{{\displaystyle \frac{\partial \mathit{\mu }}{\partial x}}}_{\mathrm{l}}-{\displaystyle \frac{D_{\mathrm{r}}{C}_{\mathrm{R}}}{{\mathit{\delta }}_x}}{{\displaystyle \frac{\partial \mathit{\varphi }}{\partial \mathit{\mu }}}}_{\mathrm{R}}{{\displaystyle \frac{\partial \mathit{\mu }}{\partial x}}}_{\mathrm{r}}.\end{array}$$

Inserting the following linearizations into Eq. (A6)

$$\begin{array}{}\text{(A7a)}& {{\displaystyle \frac{\partial \mathit{\varphi }}{\partial \mathit{\mu }}}}_{\mathrm{L}}={{\displaystyle \frac{\partial \mathit{\varphi }}{\partial \mathit{\mu }}}}_{\mathrm{C}}-{{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\varphi }}{\partial \mathit{\mu }\partial x}}}_{\mathrm{C}}{\mathit{\delta }}_x\text{(A7b)}& {{\displaystyle \frac{\partial \mathit{\varphi }}{\partial \mathit{\mu }}}}_{\mathrm{R}}={{\displaystyle \frac{\partial \mathit{\varphi }}{\partial \mathit{\mu }}}}_{\mathrm{C}}+{{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\varphi }}{\partial \mathit{\mu }\partial x}}}_{\mathrm{C}}{\mathit{\delta }}_x\end{array}$$

gives

$$\begin{array}{}\text{(A8)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{1}}{V}}{〈\sum {\mathit{\lambda }}_k{\mathit{\varphi }}_{j_k}{j}_k〉}^{**}={\displaystyle \frac{\partial \mathit{\varphi }}{\partial \mathit{\mu }}}\left({\displaystyle \frac{D_{\mathrm{l}}{C}_{\mathrm{L}}}{{\mathit{\delta }}_x}}{{\displaystyle \frac{\partial \mathit{\mu }}{\partial x}}}_{\mathrm{l}}-{\displaystyle \frac{D_{\mathrm{r}}{C}_{\mathrm{R}}}{{\mathit{\delta }}_x}}{{\displaystyle \frac{\partial \mathit{\mu }}{\partial x}}}_{\mathrm{r}}\right)\\ & -\mathrm{2}DC{\displaystyle \frac{\partial }{\partial x}}\left({\displaystyle \frac{\partial \mathit{\varphi }}{\partial \mathit{\mu }}}\right){\displaystyle \frac{\partial \mathit{\mu }}{\partial x}}.\end{array}\end{array}$$

The concentration gradient is also linearized

$$\begin{array}{}\text{(A9a)}& {\displaystyle }{C}_{\mathrm{L}}=C_{\mathrm{l}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{\partial C}{\partial x}}{\mathit{\delta }}_x,\text{(A9b)}& {\displaystyle }{C}_{\mathrm{R}}=C_{\mathrm{r}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{\partial C}{\partial x}}{\mathit{\delta }}_x.\end{array}$$

When these expressions are inserted and the limit of δ_x to zero is taken, the following partial differential equation is obtained:

$$\begin{array}{}\text{(A10)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{1}}{V}}{〈\sum {\mathit{\lambda }}_k{\mathit{\varphi }}_{j_k}{j}_k〉}^{**}=\\ & -{\displaystyle \frac{\partial \mathit{\varphi }}{\partial \mathit{\mu }}}\left[{\displaystyle \frac{\partial }{\partial x}}\left(DC{\displaystyle \frac{\partial \mathit{\mu }}{\partial x}}\right)+D{\displaystyle \frac{\partial C}{\partial x}}{\displaystyle \frac{\partial \mathit{\mu }}{\partial x}}\right]\\ & -\mathrm{2}DC{\displaystyle \frac{\partial }{\partial x}}\left({\displaystyle \frac{\partial \mathit{\varphi }}{\partial \mathit{\mu }}}\right){\displaystyle \frac{\partial \mathit{\mu }}{\partial x}}.\end{array}\end{array}$$

In the remaining steps, it will be shown that the second term on the left-hand side of Eq. (24) will cancel out with the first term on the right-hand side of Eq. (A10). Applying the product rule to $\partial \left(C\mathit{\mu }\right)/\partial t$ and substituting $\partial C/\partial t$ and $\partial \left(C\mathit{\mu }\right)/\partial t$ with Eqs. (4) and (7) results in

$$\begin{array}{}\text{(A11)}& C{\displaystyle \frac{\partial \mathit{\mu }}{\partial t}}={\displaystyle \frac{\partial }{\partial x}}\left(D{\displaystyle \frac{\partial \left(C\mathit{\mu }\right)}{\partial x}}\right)-\mathit{\mu }{\displaystyle \frac{\partial }{\partial x}}\left(D{\displaystyle \frac{\partial C}{\partial x}}\right).\end{array}$$

The product rule applied twice also implies the following result:

$$\begin{array}{}\text{(A12)}& \begin{array}{rl}{\displaystyle \frac{\partial }{\partial x}}\left(D{\displaystyle \frac{\partial \left(C\mathit{\mu }\right)}{\partial x}}\right)& ={\displaystyle \frac{\partial }{\partial x}}\left(DC{\displaystyle \frac{\partial \mathit{\mu }}{\partial x}}\right)\\ & +\mathit{\mu }{\displaystyle \frac{\partial }{\partial x}}\left(D{\displaystyle \frac{\partial C}{\partial x}}\right)+D{\displaystyle \frac{\partial C}{\partial x}}{\displaystyle \frac{\partial \mathit{\mu }}{\partial x}}.\end{array}\end{array}$$

This means that Eq. (A11) can be recast into

$$\begin{array}{}\text{(A13)}& C{\displaystyle \frac{\partial \mathit{\mu }}{\partial t}}={\displaystyle \frac{\partial }{\partial x}}\left(DC{\displaystyle \frac{\partial \mathit{\mu }}{\partial x}}\right)+D{\displaystyle \frac{\partial \mathit{\mu }}{\partial x}}{\displaystyle \frac{\partial C}{\partial x}},\end{array}$$

which matches the part between square brackets in Eq. (A10). These terms will cancel each other out when the last equation is inserted on the left-hand side and Eqs. (A3) and (A10) on the right-hand side of Eq. (24), leaving

$$\begin{array}{}\text{(A14)}& {\displaystyle \frac{\partial \left(C\mathit{\varphi }\right)}{\partial t}}={\displaystyle \frac{\partial }{\partial x}}\left(D{\displaystyle \frac{\partial \left(C\mathit{\varphi }\right)}{\partial x}}\right)-\mathrm{2}DC{\displaystyle \frac{\partial }{\partial x}}\left({\displaystyle \frac{\partial \mathit{\varphi }}{\partial \mathit{\mu }}}\right){\displaystyle \frac{\partial \mathit{\mu }}{\partial x}}.\end{array}$$

Finally, by substituting

$$\begin{array}{}\text{(A15)}& {\displaystyle \frac{\partial {\mathit{\varphi }}_q}{\partial \mathit{\mu }}}=-q{\mathit{\varphi }}_{q-\mathrm{1}},\end{array}$$

the final result shown in Table 2 is derived.