the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
A novel Eulerian model based on central moments to simulate age and reactivity continua interacting with mixing processes
Jurjen Rooze
Heewon Jung
Hagen Radtke
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 noncentral 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.
 Article
(1290 KB)  Fulltext XML

Supplement
(524 KB)  BibTeX
 EndNote
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 momentbased 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 multiG 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 spacefortime 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 centralmomentbased method is tested for age distributions and organic matter mineralization in bioturbated sediment (Sect. 4) to explore whether a momentsbased approach can effectively describe age or reactivity continua in reaction–transport models.
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
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. Lowercase and uppercase subscripts indicate evaluations at the boundaries and centers of cells, respectively. After dividing by volume V, defining $C=n/V$, and also multiplying the righthand side by ${\mathit{\delta}}_{x}^{\mathrm{2}}/{\mathit{\delta}}_{x}^{\mathrm{2}}$,
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
Applying the divergence theorem yields
which is the diffusion equation.
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_{1} and j_{4} be the fluxes that transport particles into the control volume from left and right, respectively. Similarly, let j_{2} and j_{3} 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
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_{1}=0.5f_{l}n_{L}, j_{2}=0.5f_{l}n_{C}, j_{3}=0.5f_{r}n_{C}, and j_{4}=0.5f_{r}n_{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
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,
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
The zeroth and first centralized moments are always one and zero, respectively. The variance (σ^{2}), 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 noncentral 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 qpowered differences in the control volume can be described by
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 lefthand side of Eq. 9) to the new mean age:
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 lefthand side of Eq. (10) and the first on the righthand 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
is obtained, whereby ${\mathit{\lambda}}_{k}=\pm \mathrm{1}$ depending on the direction of the flux.
2.3.1 Derivation of a partial differential equation for variance
The derivatives of variance in the Taylor series are
The only nonzero 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
or, under the assumption that $\partial \mathit{\mu}/\partial t$ is finite,
in differential form.
In the next step, the unknown fluxes on the righthand 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 higherorder 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 qpowered 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
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:
and substituting these into Eq. (15) yields
Ignoring the terms with derivatives obtained from the Taylor series for the moment, a part of the equation can be isolated
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
The remaining terms not accounted for yet are
which can also be written as
yielding in the limit of Δx→0
Therefore,
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 higherorder moments
For a finite $\partial \mathit{\mu}/\partial t$, dividing Eq. (11) by an infinitesimally small time step will drop the higherorder terms in the Taylor series, leaving
The presence of a nonzero firstorder derivative makes the derivation of the PDEs for higherorder 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.
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:
which can be obtained by applying the binomial theorem to Eq. (8). Next, we differentiate using the product rule to obtain
whereby ${\mathit{\mu}}_{x}={C}^{\mathrm{1}}\int {\mathit{\chi}}^{x}{C}_{\mathit{\chi}}\phantom{\rule{0.125em}{0ex}}\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, μ=μ_{1} 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
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 (ϕ):
The second term is further worked out by substituting
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}}_{qp}\right)/\mathrm{d}t{}_{R}$ and $\mathrm{d}\mathit{\mu}/\mathrm{d}t{}_{R}$ (Eq. 29). The first one is obtained by integrating
Following the product rule, this derivative also allows
to be solved.
3.2 Implementation of reaction kinetics and aging
Representing the production of new material with a zerothorder 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(\right)}_{\mathrm{d}\left[C{\mathit{\mu}}_{q}\right]/\mathrm{d}t}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(\right)}_{\mathrm{d}\left(C{\mathit{\mu}}_{q}\right)/\mathrm{d}t}R=\mathrm{0}$ and $\mathrm{d}\mathit{\mu}/\mathrm{d}t{}_{R}=\mathit{\mu}P{C}^{\mathrm{1}}$ into Eq. (26).
The formulation
describes firstorder 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 firstorder 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 firstorder 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.
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
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.
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 particlebased 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 finitedifferences 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^{−1}, 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.
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.
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 rootfinding procedure (Soetaert, 2009) to fit the parameter vector w.
In the example shown in Fig. 2, transport involves bioturbation and advection. The burial velocity was set to 1 mm yr^{−1}. The bioturbation coefficient had a maximum value of 10^{−10} m^{2} s^{−1} at the sediment–water interface and decreased exponentially over depth with an efolding distance set to 2 cm. The uniform distribution implemented as Dirichlet upperboundary conditions was defined by the moments C=1, $\mathit{\mu}=\mathrm{2.5}\times {\mathrm{10}}^{\mathrm{2}}$ yr^{−1}, and $\mathit{\sigma}=\mathrm{1.44}\times {\mathrm{10}}^{\mathrm{2}}$ yr^{−1}, which were also used as initial conditions throughout the domain. A nogradient condition was set as lowerboundary 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 upperboundary 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 30G 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 agedependent reactivity (Eq. 32) is specified as
whereby α prevents division by zero and may also be used as a fitting parameter. It resembles the longestablished 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 (σ^{2}=ϕ_{2}), and skewness (S=ϕ_{3}). For this, threeparameter 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:
whereby w_{1} is the multiplier, w_{2} and w_{3} denote the lower and upper limit of the distribution (outside this interval, the function evaluates to zero), and w_{4} corresponds to the mode. Given the closedform 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
which allows the distribution parameters to be calculated as follows: w_{1}=C, ${w}_{\mathrm{2}}=(\mathrm{3}\mathit{\mu}bc)/\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
whereby w_{1} serves as a multiplier to adjust the concentration, w_{2} is the shape parameter, w_{3} is the scaling parameter, and w_{4} 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 finitevolume 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 upperboundary condition is added as a Dirichlet boundary condition with a prescribed distribution. The initial values are set to the values of the upperboundary condition. A zerogradient 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).
For the simulations shown in Fig. 3, the burial velocity was set to 2 mm yr^{−1}. The bioturbation coefficient had a maximum value of 10^{−11} m^{2} s^{−1} at the top. It decreased exponentially with depth, having an efolding distance set to 3 cm, implying that the diffusivity at the bottom of the domain is effectively zero and that the zerogradient condition has a negligible effect on the results. The distributions set as upperboundary and initial conditions had moments set to C=1, μ=15 years, σ^{2}=13 yr^{2}, and S=0 yr^{3}.
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 momentbased 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 upperboundary 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.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 multiG 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 multiG 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 multiG 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 twoparameter 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 higherorder 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 finitedifferences 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 rootsolving 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, precalculated 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 momentsbased 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 noncentral 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 noncentral moments necessitate additional terms to account for aging (Delhez and Deleersnijder, 2002). While momentsbased 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 noncentral 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 multiG 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. Momentsbased simulations offer an elegant and efficient solution, while the versatile applicability of PDEs makes their implementation in various models more convenient.
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 distributionchoicerelated 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.
Continuing the derivation of the PDE for the diffusion of higher moments from Eq. (24), one can write the following equation:
where the same assumptions are made with regard to the fluxes as in the derivation for the PDE of the variance. The fluxes j_{1} and j_{4}, transporting material into the control volume, can be written as functions of the mean age at the source location,
Inserting these equations into Eq. (A1), the part not accounting for the derivatives of the Taylor series is isolated as follows:
whereby $\partial \left(C\mathit{\varphi}\right)/\partial x$ has been linearized. The terms for the derivatives can be written as follows:
Substituting this and $f/\mathrm{2}=D/{\mathit{\delta}}_{x}^{\mathrm{2}}$ into the last equation yields
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 secondorder Taylor series terms will drop. Linearizing $\partial \mathit{\mu}/\partial x$ yields
Inserting the following linearizations into Eq. (A6)
gives
The concentration gradient is also linearized
When these expressions are inserted and the limit of δ_{x} to zero is taken, the following partial differential equation is obtained:
In the remaining steps, it will be shown that the second term on the lefthand side of Eq. (24) will cancel out with the first term on the righthand 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
The product rule applied twice also implies the following result:
This means that Eq. (A11) can be recast into
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 lefthand side and Eqs. (A3) and (A10) on the righthand side of Eq. (24), leaving
Finally, by substituting
the final result shown in Table 2 is derived.
The application scripts can be accessed at https://doi.org/10.5281/zenodo.10202341 (Rooze, 2023). For any potential updates or bug fixes, please check https://git.iowarnemuende.de/rooze/DiffussionCentralMoments/ (last access: 5 December 2023). No data sets were used in this article.
The supplement related to this article is available online at: https://doi.org/10.5194/gmd1671072023supplement.
The study was conceived and designed by JR in collaboration with HJ and HR. All authors discussed the results and commented on the manuscript.
The contact author has declared that none of the authors has any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
We thank the anonymous reviewers for their constructive criticism and insightful comments.
This work was conducted within the DAM pilot mission “MGFOstsee” (grant no. 03F0848A) funded by the German Federal Ministry of Education and Research. Heewon Jung was supported by the National Research Foundation of Korea (NRF) grant funded by the South Korean government (MSIT) (grant no. 2022R1C1C1004512).
The publication of this article was funded by the Open Access Fund of the Leibniz Association.
This paper was edited by Sandra Arndt and reviewed by two anonymous referees.
Arbel, J., Lijoi, A., and Nipoti, B.: Full Bayesian inference with hazard mixture models, Comput. Stat. Data An., 93, 359–372, https://doi.org/10.1016/j.csda.2014.12.003, 2016. a
Arndt, S., Jørgensen, B. B., LaRowe, D. E., Middelburg, J., Pancost, R., and Regnier, P.: Quantifying the degradation of organic matter in marine sediments: A review and synthesis, EarthSci. Rev., 123, 53–86, https://doi.org/10.1016/j.earscirev.2013.02.008, 2013. a, b
Boudreau, B. P.: Is burial velocity a master parameter for bioturbation?, Geochim. Cosmochim. Ac., 58, 1243–1249, https://doi.org/10.1016/00167037(94)903786, 1994. a
Boudreau, B. P.: Diagenetic models and their implementation, vol. 505, Springer Berlin, https://doi.org/10.1007/9783642604218, 1997. a
Boudreau, B. P. and Ruddick, B. R.: On a reactive continuum representation of organic matter diagenesis, Am. J. Sci., 291, 507–538, https://doi.org/10.2475/ajs.291.5.507, 1991. a, b, c
Brown, P. N., Byrne, G. D., and Hindmarsh, A. C.: VODE: A variablecoefficient ODE solver, SIAM J. Sci. Stat. Comput., 10, 1038–1051, https://doi.org/10.1137/0910062, 1989. a, b
Burdige, D. J.: Preservation of organic matter in marine sediments: controls, mechanisms, and an imbalance in sediment organic carbon budgets?, Chem. Rev., 107, 467–485, https://doi.org/10.1021/cr050347q, 2007. a
Crank, J.: The mathematics of diffusion, Clarendon Press, Oxford, ISBN 13 9780198533078, 1956. a
De Borger, E., Tiano, J., Braeckman, U., Rijnsdorp, A. D., and Soetaert, K.: Impact of bottom trawling on sediment biogeochemistry: a modelling approach, Biogeosciences, 18, 2539–2557, https://doi.org/10.5194/bg1825392021, 2021. a
De Leeuw, J. W. and Largeau, C.: A Review of Macromolecular Organic Compounds That Comprise Living Organisms and Their Role in Kerogen, Coal, and Petroleum Formation, in: Organic Geochemistry: Principles and Applications, edited by: Engel, M. H. and Macko, S. A., Springer US, 23–72, https://doi.org/10.1007/9781461528906_2, 1993. a
Deleersnijder, E., Campin, J.M., and Delhez, E. J.: The concept of age in marine modelling: I. Theory and preliminary model results, J. Marine Syst., 28, 229–267, https://doi.org/10.1016/S09247963(01)000264, 2001. a, b
Delhez, É. J. and Deleersnijder, É.: The concept of age in marine modelling: II. Concentration distribution function in the English Channel and the North Sea, J. Marine Syst., 31, 279–297, https://doi.org/10.1016/S09247963(01)000665, 2002. a, b, c, d
Fick, A.: Über diffusion, Ann. Phys., 170, 59–86, https://doi.org/10.1002/andp.18551700105, 1855. a
Forbes, C., Evans, M., Hastings, N., and Peacock, B.: Statistical Distributions, John Wiley & Sons, Inc., https://doi.org/10.1002/9780470627242, 2010. a, b
Fourier, J. B. J.: Théorie analytique de la chaleur, vol. 504, Didot Paris, https://gallica.bnf.fr/ark:/12148/bpt6k1045508v (last access: 5 December 2023), 1822. a
Freitas, F. S., Pika, P. A., Kasten, S., Jørgensen, B. B., Rassmann, J., Rabouille, C., Thomas, S., Sass, H., Pancost, R. D., and Arndt, S.: New insights into largescale trends of apparent organic matter reactivity in marine sediments and patterns of benthic carbon transformation, Biogeosciences, 18, 4651–4679, https://doi.org/10.5194/bg1846512021, 2021. a, b
Gilbert, P. and Varadhan, R.: numDeriv: Accurate Numerical Derivatives, https://CRAN.Rproject.org/package=numDeriv (last access: 14 October 2023), 2016. a
Goldberg, E. D. and Koide, M.: Geochronological studies of deep sea sediments by the ionium/thorium method, Geochim. Cosmochim. Ac., 26, 417–450, https://doi.org/10.1016/00167037(62)901126, 1962. a
Guinasso Jr., N. and Schink, D.: Quantitative estimates of biological mixing rates in abyssal sediments, J. Geophys. Res., 80, 3032–3043, https://doi.org/10.1029/JC080i021p03032, 1975. a
Hasselman, B.: nleqslv: Solve Systems of Nonlinear Equations, https://CRAN.Rproject.org/package=nleqslv (last access: 14 October 2023), 2023. a
John, V., Angelov, I., Öncül, A., and Thévenin, D.: Techniques for the reconstruction of a distribution from a finite number of its moments, Chem. Eng. Sci., 62, 2890–2904, https://doi.org/10.1016/j.ces.2007.02.041, 2007. a
Jørgensen, B. B.: A comparison of methods for the quantification of bacterial sulfate reduction in coastal marine sediments: I. Measurement with radiotracer techniques, Geomicrobiol. J., 1, 11–27, https://doi.org/10.1080/01490457809377721, 1978. a, b, c
Klingbeil, K., MohammadiAragh, M., Gräwe, U., and Burchard, H.: Quantification of spurious dissipation and mixing – Discrete variance decay in a FiniteVolume framework, Ocean Model., 81, 49–64, https://doi.org/10.1016/j.ocemod.2014.06.001, 2014. a
Kuderer, M. J.: How bioturbators perturb the paleo record: From Eulerian to Lagrangian and back, Ph.D. thesis, Utrecht University, https://doi.org/10.33540/1544, 2022. a, b
Meile, C. and Van Cappellen, P.: Particle age distributions and O_{2} exposure times: timescales in bioturbated sediments, Global Biogeochem. Cy., 19, GB3013, https://doi.org/10.1029/2004GB002371, 2005. a
Meysman, F. J., Boudreau, B. P., and Middelburg, J. J.: Modeling reactive transport in sediments subject to bioturbation and compaction, Geochim. Cosmochim. Ac., 69, 3601–3617, https://doi.org/10.1016/j.gca.2005.01.004, 2005. a
Meysman, F. J., Middelburg, J. J., and Heip, C. H.: Bioturbation: a fresh look at Darwin's last idea, Trends Ecol. Evol., 21, 688–695, https://doi.org/10.1016/j.tree.2006.08.002, 2006. a
Middelburg, J. J.: A simple rate model for organic matter decomposition in marine sediments, Geochimi. Cosmochim. Ac., 53, 1577–1581, https://doi.org/10.1016/00167037(89)902391, 1989. a, b
Middelburg, J. J.: Marine carbon biogeochemistry: A primer for earth system scientists, Springer Nature, https://doi.org/10.1007/9783030108229, 2019. a
Middelburg, J. J., Soetaert, K., and Herman, P. M.: Empirical relationships for use in global diagenetic models, DeepSea Res. Pt. I:, 44, 327–344, https://doi.org/10.1016/S09670637(96)00101X, 1997. a
Narasimhan, T. N.: Fourier's heat conduction equation: History, influence, and connections, Rev. Geophys., 37, 151–172, https://doi.org/10.1029/1998RG900006, 1999. a
R Core Team: R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria, https://www.Rproject.org/ (last access: 14 October 2023), 2022. a, b
Rooze, J.: Application described in paper: A novel Eulerian model based on central moments to simulate age and reactivity continua interacting with mixing processes, Zenodo [code], https://doi.org/10.5281/zenodo.10202341, 2023. a
Rooze, J., Peterson, L., Peterson, R. N., and Meile, C.: Porewater flow patterns in surficial cold seep sediments inferred from conservative tracer profiles and early diagenetic modeling, Chem. Geol., 536, 119468, https://doi.org/10.1016/j.chemgeo.2020.119468, 2020. a
Sarmin, E. and Chudov, L.: On the stability of the numerical integration of systems of ordinary differential equations arising in the use of the straight line method, USSR Comp. Math. Math., 3, 1537–1543, https://doi.org/10.1016/00415553(63)902568, 1963. a
Soetaert, K.: rootSolve: Nonlinear root finding, equilibrium and steadystate analysis of ordinary differential equations, r package 1.6, 372 pp., ISBN 9781402086236, 2009. a
Soetaert, K. and Meysman, F.: Reactive transport in aquatic ecosystems: Rapid model prototyping in the open source software R, Environ. Modell. Softw., 32, 49–60, https://doi.org/10.1016/j.envsoft.2011.08.011, 2012. a
Soetaert, K., Petzoldt, T., and Setzer, R. W.: Solving Differential Equations in R: Package deSolve, J. Stat. Softw., 33, 1–25, https://doi.org/10.18637/jss.v033.i09, 2010. a
Tromp, T., Van Cappellen, P., and Key, R.: A global model for the early diagenesis of organic carbon and organic phosphorus in marine sediments, Geochim. Cosmochim. Ac., 59, 1259–1284, https://doi.org/10.1016/00167037(95)00042X, 1995. a
Vähätalo, A. V., Aarnos, H., and Mäntyniemi, S.: Biodegradability continuum and biodegradation kinetics of natural organic matter described by the beta distribution, Biogeochemistry, 100, 227–240, https://doi.org/10.1007/s1053301094194, 2010. a
Wang, Y. and Van Cappellen, P.: A multicomponent reactive transport model of early diagenesis: Application to redox cycling in coastal marine sediments, Geochim. Cosmochim. Ac., 60, 2993–3014, https://doi.org/10.1016/00167037(96)001408, 1996. a
Westrich, J. T. and Berner, R. A.: The role of sedimentary organic matter in bacterial sulfate reduction: The G model tested, Limnol. Oceanogr., 29, 236–249, https://doi.org/10.4319/lo.1984.29.2.0236, 1984. a
Xu, S., Liu, B., Arndt, S., Kasten, S., and Wu, Z.: Assessing globalscale organic matter reactivity patterns in marine sediments using a lognormal reactive continuum model, Biogeosciences, 20, 2251–2263, https://doi.org/10.5194/bg2022512023, 2023. a
Zindorf, M., Rooze, J., Meile, C., März, C., Jouet, G., Newton, R., Brandily, C., and Pastor, L.: The evolution of early diagenetic processes at the Mozambique margin during the last glacialinterglacial transition, Geochim. Cosmochim. Ac., 300, 79–94, https://doi.org/10.1016/j.gca.2021.02.024, 2021. a
Zonneveld, K. A. F., Versteegh, G. J. M., Kasten, S., Eglinton, T. I., Emeis, K.C., Huguet, C., Koch, B. P., de Lange, G. J., de Leeuw, J. W., Middelburg, J. J., Mollenhauer, G., Prahl, F. G., Rethemeyer, J., and Wakeham, S. G.: Selective preservation of organic matter in marine environments; processes and impact on the sedimentary record, Biogeosciences, 7, 483–511, https://doi.org/10.5194/bg74832010, 2010. a
 Abstract
 Introduction
 Derivation of partial differential equations for diffusion
 Derivation of reaction terms for partial differential equations of moments
 Applications
 Discussion
 Conclusions
 Appendix A: Derivation of diffusion PDEs for higher centralized moments
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement
 Abstract
 Introduction
 Derivation of partial differential equations for diffusion
 Derivation of reaction terms for partial differential equations of moments
 Applications
 Discussion
 Conclusions
 Appendix A: Derivation of diffusion PDEs for higher centralized moments
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement