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. 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}}$,

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.

Table 1Notation.

Download Print Version | Download XLSX

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

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

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 (σ²), 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

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:

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

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

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

Download Print Version | Download XLSX

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 }}\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

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 }}_{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

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

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.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.

Figure 1Profiles 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.

Download

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.

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 root-finding procedure (Soetaert, 2009) to fit the parameter vector w.

Figure 2A 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).

Download

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

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:

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

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

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 3Results 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.

Download

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.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.