Abstract
Reactive transport in saturated/unsaturated porous media is numerically upscaled to the spacetime scale of a hypothetical measurement through coarse grained spacetime (CGST) averages. The onedimensional reactive transport is modeled at the finegrained Darcy scale by the actual number of molecules involved in reactions which undergo advective and diffusive movements described by global random walk (GRW) simulations. The CGST averages verify identities similar to a local scale balance equation which allow us to derive expressions for the flow velocity and the intrinsic diffusion coefficient in terms of averaged microscopic quantities. The latter are further used to verify the CGSTGRW numerical approach. The upscaling approach is applied to biodegradation processes in saturated aquifers and variably saturated soils and the CGST averages are compared to classical volume averages. One finds that if the process is characterized by slow variations in time, as in homogeneous systems or in case of observations of reactive transport in heterogeneous aquifers made at large times or far away from the contaminant source, the differences between the two averages are negligible. Instead, the differences are significant if the averages are computed close to thesource at early times, in case of aquifer simulations, and can be extremely large in simulations of biodegradation in soils. In the latter case, the volume average is totally inappropriate as model for experimental measurements, leading for instance to overestimations by 100% of the CGST average.
xxx
Reactive transport in saturated/unsaturated porous media is numerically upscaled to the space–time scale of a hypothetical measurement through coarsegrained space–time (CGST) averages. The reactive transport is modeled at the finegrained Darcy scale by the actual number of molecules involved in reactions which undergo advective and diffusive movements described by global random walk (GRW) simulations. The CGST averages verify identities similar to a local balance equation which allow us to derive expressions for the flow velocity and the intrinsic diffusion coefficient in terms of averaged microscopic quantities. The latter are further used to verify the CGSTGRW numerical approach. The upscaling approach is applied to biodegradation processes in saturated aquifers and variably saturated soils and the CGST averages are compared to classical volume averages. One finds that if the process is characterized by slow variations in time, as in homogeneous reaction systems, the differences between the two averages are negligible. Instead, the differences are significant and can be extremely large in simulations of timedependent biodegradation processes in both variably saturated soils and saturated aquifers.
Authors
Nicolae Suciu
Tiberiu Popoviciu Institute of Numerical Analysis, Romanian Academy, ClujNapoca, Romania
Florin A. Radu
Department of Mathematics, University of Bergen, Allegaten 41, 5007 Bergen, Norway
Iuliu Sorin Pop
Data Science Institute, Hasselt University, Campus Diepenbeek, 3590 Diepenbeek, Belgium
Keywords
Spacetime upscaling, Global random walk, Reactive transport, Richards equation MSC: 76S05, 35K57, 86A05, 65C35
xxx
Cite this paper as:
N. Suciu, F.A. Radu, I.S. Pop, Spacetime upscaling of reactive transport in porous media, Advances in Water Resources, 176 (2023), art. no. 104443, https://doi.org/10.1016/j.advwatres.2023.104443
N. Suciu, F.A. Radu, I.S. Pop, Spacetime upscaling of reactive transport in porous media, Arxiv: https://doi.org/10.48550/arXiv.2112.10692
About this paper
Journal
Advances in Water Resources
Publisher Name
Elsevier Ltd.
Print ISSN
03091708
Online ISSN
Not available yet.
Google Scholar Profile
Spacetime upscaling of reactive transport in porous media
Abstract
Reactive transport in saturated/unsaturated porous media is numerically upscaled to the spacetime scale of a hypothetical measurement through coarsegrained spacetime (CGST) averages. The reactive transport is modeled at the finegrained Darcy scale by the actual number of molecules involved in reactions which undergo advective and diffusive movements described by global random walk (GRW) simulations. The CGST averages verify identities similar to a local balance equation which allow us to derive expressions for the flow velocity and the intrinsic diffusion coefficient in terms of averaged microscopic quantities. The latter are further used to verify the CGSTGRW numerical approach. The upscaling approach is applied to biodegradation processes in saturated aquifers and variably saturated soils and the CGST averages are compared to classical volume averages. One finds that if the process is characterized by slow variations in time, as in homogeneous reaction systems, the differences between the two averages are negligible. Instead, the differences are significant and can be extremely large in simulations of timedependent biodegradation processes in both variably saturated soils and saturated aquifers.
Keywords:
Spacetime upscaling, Global random walk, Reactive transport, Richards equation MSC: 76S05, 35K57, 86A05, 65C351 Introduction
Concentrations measured in soil experiments as well as in samples pumped from monitoring wells or collected by multilevel samplers are necessarily space averages over the support scale of the sampling procedure and measuring device [Mailloux et al.,, 2003, Teutsch et al.,, 2000, BayerRaich et al.,, 2004]. The sampling volume enters as an important parameter in analytical and numerical models because of its impact on the moments and the probability distribution of the uncertain concentration transported in heterogeneous natural porous media [Andričević,, 1998, Schwede et al.,, 2008, Tonina and Belin.,, 2008, Srzic et al.,, 2013, de Barros and Fiori,, 2014].
The spatial average over the sampling volume is accounted for in different ways depending on modeling approaches. In analytical and semianalytical approaches, the average concentration is an integral of a continuous field of concentrations [Andričević,, 1998, Schwede et al.,, 2008] or related probability densities [de Barros and Fiori,, 2014]. The same procedure is applied to both resident and fluxweighted concentrations [Schwede et al.,, 2008]. When the purpose is to estimate probability distributions, the resident concentration is preferred [Schwede et al.,, 2008, de Barros and Fiori,, 2014]. The resident concentrations solve local balance equations for reactive transport which are the starting point in deriving evolution equations for the concentration probability density function [Suciu et al.,, 2015] and the corresponding concentration moments equations [Schüler et al.,, 2016]. These equations describe the statistics of the concentration without considering a sampling volume, which corresponds to a point sampling in volume averaging approaches [Tonina and Belin.,, 2008, Srzic et al.,, 2013]. However, by averaging the solution of the local balance equation with a spatial filter, a filtered density function approach is obtained, which explicitly takes into account the sampling volume, determined by the width of the spatial filter [Suciu et al.,, 2016].
In the case of discrete descriptions of the transport process, such as sequential particle tracking [Caroni and Fiorotto,, 2005, Srzic et al.,, 2013, Tonina and Belin.,, 2008, Wright et al.,, 2021] and global random walk (GRW) simulations [Suciu et al.,, 2016, Suciu,, 2019, Suciu et al.,, 2021, Suciu and Radu,, 2022], concentrations are measured by counting particles in the sampling volume. When using moderate numbers of particles, as in particle tracking approaches, the traditional counting particles produces nonsmooth fields. This issue is overcome in kernel density estimations [FernàndezGarcia and SanchezVila,, 2011, SoleMari and FernàndezGarcia,, 2018, SoleMari et al.,, 2019] by weighting the contributions of the individual particles with smooth kernel functions. The resulting kernel estimates are smooth concentration fields and concentration derivatives/gradients. The parameter of the kernel function is a measure of the support volume which defines the region over which a particle influences the estimate. The optimal volume support can be evaluated as function of the number of particles and it has been shown that it goes to zero for infinite number of particles [FernàndezGarcia and SanchezVila,, 2011]. The need of an optimal support volume is thus related to the relatively small numbers of particles used in particle tracking approaches, which is limited by the computational resources. In GRW methods which, instead of tracking individual particles, move groups of particles over the sites of a regular lattice, particles can be associated to molecules and a mole of solute is represented by Avogrado’s number [Suciu and Radu,, 2022].
The relation between the two spatial averaging procedures, by integrating continuous fields and by counting particles, is not trivial. In fact, this is equivalent to passing from a microscopic level of description to a macroscopic one. Such connections can be rigorously established if the spatial sampling is replaced by a spacetime average [Suciu,, 2001, Vamoş,, 2007].
Considering the time scale along with the support volume is imposed by the obvious fact that any measurement is also a time average and hydrological measurements are no exception [Destouni and Graham,, 1997, Andričević,, 1998]. Even if some flow and transport conditions (e.g. small spatial scale of concentration fluctuations [Andričević,, 1998] or steady state transport [Schwede et al.,, 2008]) justify a sampling volume approach which disregards the measurement time scale, in most cases both the spatial and temporal scales specific to the observation method [BayerRaich et al.,, 2004, Berkowitz,, 2021, Destouni and Graham,, 1997] have to be considered.
The influence of the observation method on the expected value and the variance of the locally measured concentrations has been modeled in [Destouni and Graham,, 1997] by time averages of Lagrangian concentration representations. The proposed methodology can be used for either unsaturated or saturated flow conditions and for different types of observation methods. The temporal scale may be related to a spatial scale, as for instance for constant water volume sampling or for constant soil sample length. But the temporal scale may also be independent of the spatial scale, as in case of monitoring wells, where it is related to the pumping time [BayerRaich et al.,, 2004], or in lysimeter experiments [Dabrowska and Rykala,, 2021], where the relevant scale is the effluent collection time. The analysis based on temporal averages presented in [Destouni and Graham,, 1997] does not consider the effect of the spatial averaging and the results are consistent with those reported in previously published volume sampling studies.
A theoretical approach based on spatial and temporal averages was developed to transform microscale balance equations into macroscopic equations [He and Sykes,, 1996]. The basic assumption in this approach is the existence of differential equations which describe the transport of the volumetric density of an extensive physical quantity at the microscale. The averaging domain consists of the set product between a representative elementary volume which may vary in space and time and a variable timescale. The spacetime averaging domain can thus be thought of as a generalization of the classical “fluid molecule”. Assuming the existence of the solution of the microscale equation, averaging theorems which relate derivatives of averages to averages of derivatives are proved and further used to derive the macroscopic equation for multiphase transport in porous media [Cushman,, 1983, He and Sykes,, 1996]. This equation is not closed and constitutive relations are needed to formulate wellposed problems.
Most recent papers on averaging procedures for porous media do not consider the temporal averaging, focusing mainly on volume averaging theorems for multiphase systems (see [Gray and Miller.,, 2013] and references therein). Spatial averaging procedures are used, for instance, to formulate macroscopic equations for multiphase transport [Gray et al.,, 2015], to understand how the internal structure of a twofluid porous medium system influences the macroscale state [McClure et al.,, 2017], to account for dependencies between the size of the averaging volume and the macroscale in derivations of Darcy’s law from NavierStokes equations [Nordbotten et al.,, 2007], as well as in upscaling transport with Monod reactions to obtain effective equations for the macroscale [Heße et al.,, 2009].
Important achievements in deriving upscaled models for porous media were obtained by using homogenization techniques [Hornung,, 1996] which, unlike direct volume or spatiotemporal averaging, lead to closed upscaled equations with macroscopic coefficients completely determined by solving unitcell problems [Auriault and Adler,, 1995]. Examples encompass but are not limited to, derivation of the Darcy scale model for precipitation and dissolution in porous media [Kumar et al.,, 2016], modeling injection of cold water in geothermal reservoirs [Bringedal et al.,, 2016], upscaling flow by homogenization methods for nonperiodic and random parameters of the subsurface system [Nolen et al.,, 2008]. The concept of spatiotemporal upscaling in the context of homogenization by multiplescale expansions, initially developed to study linear equations with coefficients oscillating in both space and time [Bensoussan et al.,, 2011], has been extended to nonlinear diffusion equations [Akagi and Oka,, 2020] and to reactive transport in porous media driven by oscillating pore water velocity [Van Duijn and van der Zee,, 2018]. In this approach, spatial and temporal scales are necessarily related, which could be a limitation of the spatiotemporal homogenization from the point of view of engineering applications. Another inconvenience is that the spatial and temporal scales are related to the inverse of the small parameter of the homogenization problem and the macroscopic description obtained when this parameter goes to zero cannot be directly related to the scales of measurement instruments and methods [Destouni and Graham,, 1997, BayerRaich et al.,, 2004] or to the spatial and temporal scales of the hydrological processes [Skøien et al.,, 2003].
In both direct spatiotemporal averaging and homogenization approaches the “microscopic” state is given by continuous fields. The alternative “counting particles” averaging approach requires a microscopic description given by trajectories of physical or computational particles. Assuming that the heterogeneity of the subsurface system can be modeled by random functions, a “microscopic” description corresponding to the local, Darcy scale of the transport process is given by trajectories of diffusion in random velocity fields [Suciu,, 2019] modeled numerically by particle tracking [Wright et al.,, 2021] or GRW algorithms [Suciu et al.,, 2021]. A general method to obtain a macroscopic description at spatiotemporal scales that can be directly related to those of the measurements and hydrological observations is the coarsegrained (CG) description of corpuscular systems through spacetime (ST) averages [Vamoş et al.,, 1996]. The CGST approach also yields unclosed macroscopic equations but, since a complete microscopic description is available, it provides a frame to investigate the structure of the constitutive relations through numerical simulations [Vamoş et al.,, 1997].
Upscaling reactive transport in subsurface hydrological systems aims at obtaining models that can account for small scale processes in predicting the fate of contaminants at larger scales and in designing remediation strategies for contaminated sites. Solutions proposed so far, based on either Eulerian approaches [Porta et al.,, 2012] or on particle tracking transport models combined with Eulerian descriptions of the reactions [Wright et al.,, 2021] are essentially volume averaging methods. However, it has been argued that temporal fluctuations of reacting species concentrations due to the extrinsic noise resulted from the variability of the chemical environment and of the transport process [Aquino and Dentz,, 2017] propagate from local to larger spatial scales and at least an implicit temporal averaging must also be included in upscaling reactive transport [Berkowitz,, 2021]. It is further argued in [Berkowitz,, 2021] that the temporal information can be included through a continuoustime random walk approach, with a powerlaw distribution of transition times calibrated from measurements, which accounts for temporal effects over a range of time scales.
The CGST method offers the perspective of explicitly relating the upscaled concentrations to the spatiotemporal scales of interest. The finegrained concentrations are represented microscopically by the actual number of molecules involved in reactions, their advectivediffusive movement being described by simple discretetime random walk processes modeled with GRW algorithms [Suciu and Radu,, 2022]. This molecular description implicitly accounts for spatiotemporal fluctuations produced by inhomogeneities of reaction and transport parameters. Moreover, by using integrated flow and transport GRW solvers for unsaturated/saturated porous media [Suciu et al.,, 2021] the CGST approach can be applied to upscaling contaminant transport in both soil and aquifer systems.
In this article we propose a CGST upscaling approach using microscopic descriptions consisting of one and twodimensional GRW simulations of the reactive transport at Darcy scale. General statistical mechanics results on CGST averaging are presented in Section 2. In Section 3, theoretical formulas for the flow velocity and the intrinsic diffusion coefficient are derived as corollaries of the general result and are used to verify the Matlab codes designed for the numerical spacetime upscaling. Further, in Section 4 the CGST approach is applied to various onedimensional problems of reactive transport solved by GRW algorithms, with emphasis on biodegradation processes in saturated and unsaturated natural porous media. Twodimensional codes for CGST upscaling are verified and applied to biodegradation in variably saturated soils in Section 5. The main results and summarized in Section 6. The proof of main result presented in Section 2 is given in A and B illustrates the influence of the spacetime scale on the discrepancy between volume and CGST averages. The codes implemented in Matlab used in this study are stored in the Git repository SpaceTimeUpscaling [Suciu et al., 2023a, ].
2 Coarsegrained spacetime averages
The CGST method allows continuous modeling of corpuscular systems, without the need to use the vague notion of “fluid particle” which is small enough to be infinitesimal with respect to the size of the macroscopic system, but large enough to contain many molecules in local thermodynamic equilibrium. Continuous macroscopic fields and balance equations can be derived under the only assumption that a system of particles (e.g. molecules and other physical particles, or computational particles used in GRW and particle tracking simulations) possesses a kinematic microscopic description through piecewise analytic time functions. A general proof of this result has been obtained in [Vamoş et al.,, 1996] by considering CGST averages over spheres and symmetric time intervals. Proofs for onedimensional systems were presented in [Vamoş et al.,, 1997, 2000] and a general result for averaging on hypercubes in the phase space of the statistical mechanics, for conserved number of phase points, was achieved in [Vamoş,, 2007]. In this article we consider a general kinematic description for systems of particles that can be created and destroyed, with piecewise analytical trajectories and discontinuous velocities at a finite number of points. Under these conditions, we prove in A the following properties of the CGST averages over $d$dimensional cubes and symmetric time intervals.
If a discrete system of $\mathcal{N}$ particles is described by piecewise analytic time functions ${\phi}_{i}(t):[0,T]\u27fc\mathbb{R}$, $i=1,\mathrm{\dots},\mathcal{N}$, then
1) there exists a macroscopic description of the same system given by almost everywhere (a.e.) continuous fields through coarsegrained spacetime averages
$$\u27e8\phi \u27e9(\mathbf{x},t;a,\tau )=\frac{1}{2\tau {(2a)}^{d}}\sum _{i=1}^{\mathcal{N}}{\int}_{t\tau}^{t+\tau}{\phi}_{i}({t}^{\prime})\prod _{\alpha =1}^{d}{\chi}_{\alpha}({x}_{\alpha i}({t}^{\prime}))d{t}^{\prime},$$  (1) 
where ${x}_{\alpha i},\alpha =1,\mathrm{\dots},d$, are the components of the $i$th particle trajectory, $$ is the half length of the temporal averaging interval, ${(2a)}^{d}$ is the volume of the open $d$dimensional cube $$, and ${\chi}_{\alpha}$ is the characteristic function of the interval $({x}_{\alpha}a,{x}_{\alpha}+a)$,
2) $\u27e8\phi \u27e9(\mathbf{x},t;a,\tau )$ has continuous partial derivatives a.e. in ${\mathbb{R}}^{d}\times (\tau ,T\tau )$ and satisfies the identity
$${\partial}_{t}\u27e8\phi \u27e9+\nabla \cdot \u27e8\phi \bm{\xi}\u27e9=\u27e8\frac{d\phi}{dt}\u27e9+\delta \phi ,$$  (2) 
where ${\xi}_{\alpha i}=d{x}_{\alpha i}/dt$ are the velocity components of the $i$th particle and $\delta \phi $ accounts for discontinuous variations produced when a particle is created or consumed in chemical reactions or by finite jumps of the function ${\phi}_{i}(t)$.
For ${\phi}_{i}\equiv 1$, $i=1\mathrm{\dots}\mathcal{N}$, the CGST average (1) counts the number of particles inside the cube $C(\mathbf{x},a)$ during a time interval of length $2\tau $, that is, $\u27e81\u27e9$ is an a.e. continuous approximation of the macroscopic concentration field. The average over the number of particles defined by $\overline{\bm{\xi}}=\u27e8\bm{\xi}\u27e9/\u27e81\u27e9$, if $\u27e81\u27e9>0$, and by $\overline{\bm{\xi}}=0$, if $\u27e81\u27e9=0$ provides an a.e. continuous approximation of the macroscopic velocity field [Vamoş et al.,, 1996, 2000]. In particular, if the particles move with the constant velocity $\mathbf{u}$, according to (1), $\overline{\bm{\xi}}=\mathbf{u}$. In the absence of chemical reactions and jump discontinuities the term $\delta \phi $ in (2) vanishes and one obtains a continuity equation,
$${\partial}_{t}\u27e81\u27e9+\nabla \cdot (\u27e81\u27e9\overline{\bm{\xi}})=0.$$  (3) 
For $\tau \u27f60$ the CGST average (1) with ${\phi}_{i}\equiv 1$ gives the estimation of the concentration by the usual volume averaging formula. The essential advantage of (1) is that the additional time average yields a.e. continuous fields obeying the identity (2), which provides a frame to close the balance equations [Vamoş et al.,, 1997, 2000].
In the framework of the statistical mechanical theory of transport processes [Irving and Kirkwood,, 1950], the continuous concentration field $c(\mathbf{x},t)$ at the scale of the mathematical continuum, hereafter called “finegrained concentration”, is defined as an average $M$ over the ensemble of realizations of the transport process,
$$c(\mathbf{x},t):=M\left[\sum _{i=1}^{\mathcal{N}}\delta ({\mathbf{X}}_{i}(t,\omega )\mathbf{x})\right],$$ 
where $\delta $ is the Dirac function and ${\mathbf{X}}_{i}(t,\omega )$ denotes the realization $\omega $ of the random trajectory. Applying the averaging operator $M$ to (1), for ${\phi}_{i}\equiv 1$, with ${x}_{\alpha i}(t)={X}_{\alpha i}(t,\omega )$, and expressing the characteristic function of the cube in (1) with the aid of the Dirac functional
$$\prod _{\alpha =1}^{d}{\chi}_{\alpha}({X}_{\alpha i}(t))={\chi}_{C(\mathbf{x},a)}({\mathbf{X}}_{i},t)={\int}_{{\mathbb{R}}^{d}}{\chi}_{C(\mathbf{x},a)}({x}^{\prime})\delta ({\mathbf{X}}_{i}(t,\omega ){\mathbf{x}}^{\prime})\mathit{d}{\mathbf{x}}^{\prime},$$ 
one obtains the following result.
The ensemble average of the coarsegrained average $\u27e81\u27e9(\mathbf{x},t;a,\tau )$ is an average of the finegrained concentration $c(\mathbf{x},t)$ over the cube $C(\mathbf{x},a)$ and the time interval $(t\tau ,t+\tau )$ given by
$$M[\u27e81\u27e9](\mathbf{x},t;a,\tau )=\frac{1}{2\tau {(2a)}^{d}}\underset{t\tau}{\overset{t+\tau}{\int}}\mathit{d}{t}^{\prime}\underset{C(\mathbf{x},a)}{\int}c({\mathbf{x}}^{\prime},{t}^{\prime})\mathit{d}{\mathbf{x}}^{\prime}.$$  (4) 
The ensemble average $M[\u27e81\u27e9](\mathbf{x},t;a,\tau )$ defines the continuous concentration field $c(\mathbf{x},t;a,\tau )$ as observed at the spatiotemporal scale $(a,\tau )$. Equation (4) is a particular case of the general relation between the smooth continuous field and the CGST average of its microscopic description [Suciu,, 2001, Eq. (4.6)]. For $\tau \u27f60$, Eq. (4) gives the relation between the “integrating continuous fields” and “counting particles” averaging procedures in sampling volume approaches.
The use of CGST averages (1) in numerical simulations is conditioned by the piecewise analyticity of $\phi (t)$. Random walks on lattices are eligible because the corresponding trajectories ${\mathbf{x}}_{i}(t)$ and velocities ${\bm{\xi}}_{i}$ are piecewise analytic functions with numerable discontinuity points at the jump instants [Vamoş et al.,, 2000]. In GRW simulations, the sum over the number of particles in (1) is equivalent to summing up the contributions of the groups of particles jumping between the same lattice sites in the interior of the cube $C(\mathbf{x},a)$, in the time interval $(t\tau ,t+\tau )$.
3 Verification of the CGST averaging procedure
The goal of the spacetime upscaling will be achieved by performing CGST averages on microscopic descriptions consisting of simulations of transport in porous media performed with GRW algorithms. In this section, we consider the onedimensional case and we verify the reliability of the averaging procedure applied to simulations of diffusion with both the unbiased GRW and its biased version, BGRW.
3.1 Implications of the CGST upscaling
In the onedimensional case, the CGST average (1) becomes
$$\u27e8\phi \u27e9(x,t;a,\tau )=\frac{1}{4\tau a}\sum _{i=1}^{\mathcal{N}}\underset{t\tau}{\overset{t+\tau}{\int}}{\phi}_{i}({t}^{\prime})\chi ({x}_{i}({t}^{\prime}))\mathit{d}{t}^{\prime}$$  (5) 
and the continuity equation (3) takes the form
$${\partial}_{t}\u27e81\u27e9+{\partial}_{x}(\u27e81\u27e9\overline{\xi})=0.$$  (6) 
Considering ${\phi}_{i}={x}_{i}$ and using (2), the advective flux in (6) can be expressed as
$$\u27e81\u27e9\overline{\xi}=\u27e8\xi \u27e9=\u27e8\frac{dx}{dt}\u27e9={\partial}_{t}\u27e8x\u27e9+{\partial}_{x}\u27e8x\xi \u27e9.$$ 
Further, by assuming $\u27e81\u27e9>0$ and considering the average $\overline{x}=\u27e8x\u27e9/\u27e81\u27e9$ of particle positions in the averaging domain $(xa,x+a)\times (t\tau ,t+\tau )$ and the continuity equation (6), we get
$${\partial}_{t}\u27e8x\u27e9={\partial}_{t}(\u27e81\u27e9\overline{x})=\u27e81\u27e9{\partial}_{t}\overline{x}+\u27e81\u27e9\overline{\xi}{\partial}_{x}\overline{x}{\partial}_{x}(\u27e81\u27e9\overline{x}\overline{\xi})$$ 
and the continuity equation can be reformulated as a FokkerPlanck equation,
$${\partial}_{t}\u27e81\u27e9+{\partial}_{x}[\u27e81\u27e9({\partial}_{t}\overline{x}+\overline{\xi}{\partial}_{x}\overline{x})]={\partial}_{x}^{2}[(\overline{x}\overline{\xi}\overline{x\xi})\u27e81\u27e9].$$  (7) 
The positivity of the coefficient $(\overline{x}\overline{\xi}\overline{x\xi})$ in (7) is a “microscopic criterion of irreversibility” for the thermodynamic system of the $\mathcal{N}$ particles [Suciu,, 2001]. However, the drift coefficient $({\partial}_{t}\overline{x}+u{\partial}_{x}\overline{x})$ differs from the Eulerian drift $u$ and Eq. (7) cannot be generally identified with the advectiondiffusion equation valid at the macroscopic scale. Instead, assuming the thermodynamic equilibrium of the system of particles, ${\partial}_{t}\u27e8x\u27e9=0$, and $u=0$, which implies ${\partial}_{t}c=0$ according to (6), Eq. (7) reduces to
$${\partial}_{x}^{2}(D\u27e81\u27e9)=0,$$  (8) 
Equation (8) is a stationary diffusion equation with diffusion coefficient $D=\overline{x\xi}=\u27e8x\xi \u27e9/\u27e81\u27e9$, which under the assumption of thermodynamic equilibrium must be constant. Hence, the reformulation of the continuity equation and the assumption of thermodynamic equilibrium highlight the selfdiffusion process with “intrinsic diffusion coefficient” $D$ of the system of particles of the microscopic description.
The relation defining the intrinsic diffusion coefficient $D=\u27e8x\xi \u27e9/\u27e81\u27e9$, as well as the average defining the macroscopic velocity $u=\u27e8\xi \u27e9/\u27e81\u27e9$ will be used in the following to verify the numerical approach for computing CGST averages over microscopic descriptions consisting of systems of random walkers.
3.2 GRW/BGRW algorithms
The GRW algorithms can be understood as weak schemes for Itô stochastic differential equations [Suciu,, 2019]. In the above onedimensional case, a weak scheme for solving the equation $d{X}_{t}=udt+\sqrt{2D}dW$, where $W$ is the standard Wiener process, can be obtained by using the discrete random walk process ${X}_{k\mathrm{\Delta}t}$ with increments ${X}_{(k+1)\mathrm{\Delta}t}{X}_{k\mathrm{\Delta}t}=(\stackrel{~}{u}+\zeta d)\mathrm{\Delta}x$ consisting of jumps over $(\stackrel{~}{u}+\zeta d)$ steps $\mathrm{\Delta}x$ on a regular lattice. The integer number $d$ is the amplitude of the diffusive displacements, $\stackrel{~}{u}=\lfloor u\mathrm{\Delta}t/\mathrm{\Delta}x\rfloor $ is the truncated advective displacement, $\zeta $ is a random variable taking three discrete values with probabilities $P(\zeta =\pm 1)=r/2$, $P(\zeta =0)=(1r)$, and the real parameter $r=2D\mathrm{\Delta}t/{(d\mathrm{\Delta}x)}^{2}$ is the dimensionless diffusion coefficient. The probability density of the Itô process ${X}_{t}$, i.e. the normalized concentration of random walkers $c(x,t)$ solving the FokkerPlanck equation ${\partial}_{t}c+{\partial}_{x}(uc)=D{\partial}_{x}^{2}c$, is approximated by counting the number of random walkers at lattice sites $l$, $c(l\mathrm{\Delta}x,k\mathrm{\Delta}t)\approx {n}_{l,k}/\mathcal{N}$, where $\mathcal{N}$ is the total number of random walkers. As in Section 2 above, by $c(x,t)$ we denote the finegrained concentration.
The unbiased GRW algorithm distributes the random walkers over lattice sites according to
$${n}_{l,k}=\delta {n}_{l+\stackrel{~}{u}l,k}+\delta {n}_{l+\stackrel{~}{u}dl,k}+\delta {n}_{l+\stackrel{~}{u}+dl,k},$$  (9) 
where the numbers $\delta n$ of random walkers are binomial random variables with statistics specified by the parameter $r$. They represent the number of particles which at time $k$ undergo jumps from a lattice site $l$ to $l+\stackrel{~}{u}$, $l+\stackrel{~}{u}d$, and $l+\stackrel{~}{u}+d$, respectively.
An alternative weak scheme for the above Itô equation uses the random walk process with increments ${X}_{(k+1)\mathrm{\Delta}t}{X}_{k\mathrm{\Delta}t}=\zeta \mathrm{\Delta}x$ and biased jump probabilities $P(\zeta =\pm 1)=(r\pm u)/2$, $P(\zeta =0)=(1r)$. This results in a biased algorithm, BGRW, which distributes the random walkers on firstneighbor lattice sites according to
$${n}_{l,k}=\delta {n}_{ll,k}+\delta {n}_{l1l,k}+\delta {n}_{l+1l,k}.$$  (10) 
The restriction $u\le r$ implies an upper bound for the local Péclet number, $\text{P\xe9}=u\mathrm{\Delta}x/D\le 2$, which imposes smaller $\mathrm{\Delta}x$ and $\mathrm{\Delta}t$, rendering the BGRW algorithm slower than the unbiased GRW algorithm.
For a general presentation of the GRW/BGRW algorithms we refer to [Suciu,, 2019, Chap. 3]. Implementation details, for more complex problems consisting of onecomponent reactive transport fully coupled to saturated/unsaturated flows governed by the degenerate Richards equation can be found in [Suciu et al.,, 2021, Sect. 4.2], and for transport with multispecies Monod reactions in [Suciu and Radu,, 2022, Sect. 3].
3.3 Averaging procedure and verification
The CGST average (5) on the GRW/BGRW simulations is computed by summing up contributions from groups $\delta n$ of random walkers as follows. Let ${N}_{a}$ be the number of lattice sites inside the open interval $(xa,x+a)$. If the conditions ${X}_{k\mathrm{\Delta}t}x\le a$ and $$ are fulfilled at every occupied site, the CGST average of the microscopic quantity $\phi $ is computed by
$$\u27e8\phi \u27e9(x,t;a,\tau )=\frac{\mathrm{\Delta}t}{4\tau a}\sum _{k=(t\tau )/\mathrm{\Delta}t}^{(t+\tau )/\mathrm{\Delta}t}\sum _{i=1}^{{N}_{a}}{\mathrm{\Phi}}_{i,k}(\phi ),$$  (11) 
where ${\mathrm{\Phi}}_{i,k}(\phi )$ are contributions of the quantity $\phi $ to the CGST average.
As microscopic quantities we consider in the following ${\phi}_{k}=1$, for the computation of the CGST concentration $\u27e81\u27e9$, ${\phi}_{k}={X}_{(k+1)\mathrm{\Delta}t}$, for the mean particle positions $\u27e8x\u27e9$, and ${\phi}_{k}=({X}_{(k+1)\mathrm{\Delta}t}{X}_{k\mathrm{\Delta}t})/\mathrm{\Delta}t$, for the mean velocity $\u27e8\xi \u27e9$. The computation of the intrinsic diffusion coefficient $D=\u27e8x\xi \u27e9/\u27e81\u27e9$ defined during the derivation of the stationary diffusion equation (8) will be also used in the following to verify the correctness of the averaging procedure. Since $D$ is a property of the thermodynamic system, it does not depend on the flow regime and can be computed as well for microscopic descriptions consisting of GRW simulations of advectiondiffusion processes. To do that, we have to consider the velocity $\xi $ of the random walk jumps. In the case of the microscopic description provided by the unbiased GRW algorithm, $\xi $ computed by disregarding the advective displacement $\stackrel{~}{u}$ is associated to the middle of the jump interval and $\u27e8x\xi \u27e9$ is computed for the microscopic quantity ${\phi}_{k}=({X}_{k\mathrm{\Delta}t}+\stackrel{~}{u}\mathrm{\Delta}x+\frac{1}{2}\zeta d\mathrm{\Delta}x)(\zeta d\mathrm{\Delta}x)/\mathrm{\Delta}t$. When using the BGRW algorithm (10), where the bias of the random walk jump probabilities accounts for advection, $\u27e8x\xi \u27e9$ is computed with binomial random variables $\delta {n}_{l\pm 1\mid l,k}$ representing the number of random walkers jumping on neighboring sites with unbiased probabilities $\pm r/2$, where $r$ is the same parameter as that used to compute biased jump probabilities $(r\pm u)/2$ in case of nonvanishing advection velocity $u$, and the corresponding microscopic quantity is given by ${\phi}_{k}=({X}_{k\mathrm{\Delta}t}+\frac{1}{2}\zeta \mathrm{\Delta}x)(\zeta \mathrm{\Delta}x/\mathrm{\Delta}t)$.
Using (9) and (10), the CGST averages considered in this article are computed according to (11) with the functions ${\mathrm{\Phi}}_{i,k}(\phi )$ listed below:

1.
the concentration $\u27e81\u27e9$
${\mathrm{\Phi}}_{i,k}(\phi )={n}_{i,k+1}$

2.
the mean particle positions $\u27e8x\u27e9$

(a)
GRW algorithm: ${\mathrm{\Phi}}_{i,k}(\phi )=[(i+\stackrel{~}{u})\mathrm{\Delta}x]{n}_{i,k}$

(b)
BGRW algorithm: ${\mathrm{\Phi}}_{i,k}(\phi )=(i\mathrm{\Delta}x)\delta {n}_{ii,k}+[(i1)\mathrm{\Delta}x]\delta {n}_{i1i,k}+[(i+1)\mathrm{\Delta}x]\delta {n}_{i+1i,k}$

(a)

3.
the mean velocity $\u27e8\xi \u27e9$

(a)
GRW algorithm: ${\mathrm{\Phi}}_{i,k}(\phi )=[\stackrel{~}{u}(\mathrm{\Delta}x/\mathrm{\Delta}t)]{n}_{i,k}$

(b)
BGRW algorithm: ${\mathrm{\Phi}}_{i,k}(\phi )=(\mathrm{\Delta}x/\mathrm{\Delta}t)\delta {n}_{i1i,k}+(\mathrm{\Delta}x/\mathrm{\Delta}t)\delta {n}_{i+1i,k}$

(a)

4.
the product positionvelocity $\u27e8x\xi \u27e9$

(a)
GRW algorithm:
${\mathrm{\Phi}}_{i,k}(\phi )=[(i+\stackrel{~}{u}\frac{1}{2}d)\mathrm{\Delta}x](d\mathrm{\Delta}x/\mathrm{\Delta}t)\delta {n}_{i+\stackrel{~}{u}di,k}+[(i+\stackrel{~}{u}+\frac{1}{2}d)\mathrm{\Delta}x](d\mathrm{\Delta}x/\mathrm{\Delta}t)\delta {n}_{i+\stackrel{~}{u}+di,k}$

(b)
BGRW algorithm:
${\mathrm{\Phi}}_{i,k}(\phi )=(i\frac{1}{2})\mathrm{\Delta}x(\mathrm{\Delta}x/\mathrm{\Delta}t)\delta {n}_{i1i,k}^{(unbiased)}+(i+\frac{1}{2})\mathrm{\Delta}x(\mathrm{\Delta}x/\mathrm{\Delta}t)\delta {n}_{i+1i,k}^{(unbiased)}$

(a)
The averaging procedure is verified by computing the flow velocity and the intrinsic diffusion coefficient with the expressions derived by CGST averaging in Section 3.1. We shall consider the fine scale transport process described by the advectiondiffusion equation with constant coefficients $D$ and $u$,
$${\partial}_{t}c+u{\partial}_{x}(c)=D{\partial}_{x}^{2}c.$$  (12) 
The finegrained concentration $c$ solving (12) is approximated numerically with biased and unbiased GRW algorithms, which also provide microscopic descriptions of the physical system. As follows from Eq. (5), in this case the macroscopic flow velocity equals the constant drift $u$. Also, one expects that the intrinsic diffusion coefficient equals the diffusion parameter $D$ of the advectiondiffusion process.
For illustration and verification purposes, we consider the onedimensional spatial domain $\mathrm{\Omega}=(0,1)$, a time interval of length $T=1$, and GRW/BGRW lattices of constant $\mathrm{\Delta}x$ and $\mathcal{L}=1/\mathrm{\Delta}x+1$ sites, denoted by $l=1,\mathrm{\dots},\mathcal{L}$. We impose noflux conditions ${n}_{1,k}={n}_{2,k}$ and ${n}_{\mathcal{L},k}={n}_{\mathcal{L}1,k}$ which ensure a steadystate solution of Eq. (12). The initial condition consists of $\mathcal{N}={10}^{24}$ random walkers uniformly distributed on the lattice, that is, a constant concentration ${n}_{l,k}/\mathcal{N}=1/\mathcal{L}$ mole per lattice site. Under these conditions, according to (11), one obtains the normalized CGST average concentration $\u27e81\u27e9/\mathcal{N}={N}_{a}/(2a\mathcal{L})$, which does not depend on $\tau $ and equals the spatial average over the interval $(xa,x+a)$. With $a=0.05$ and $\tau =0.1$, the CGST averages are computed over 9 disjoint spatial intervals of length $2a$ and three time intervals of length $2\tau $ centered at the sampling times $t=$0.1, 0.5, and 0.9. The lattice constant is progressively decreased from $\mathrm{\Delta}x=5\cdot {10}^{3}$ to $\mathrm{\Delta}x=3.33\cdot {10}^{4}$ and the corresponding time steps $\mathrm{\Delta}t$ are obtained from the condition $r=2D\mathrm{\Delta}t/{(d\mathrm{\Delta}x)}^{2}\le 1$, with $d=1$ in case of BGRW algorithm, and with both $d=1$ and $\stackrel{~}{u}=1$ in case of unbiased GRW (See Section 3.2). The constant velocity is set to $u=1$ and the diffusion coefficient is set to $D=0.0001$.
It has been found that the steady state, as indicated by a CGST concentration almost identical to that given by the spatial average, is reached for $\mathrm{\Delta}x\le 5\cdot {10}^{4}$ in case of the BGRW microscopic description (see Figs. 2 and 2) and for $\mathrm{\Delta}x\le 3.33\cdot {10}^{4}$ in case of the GRW description. However, as seen in Figs. 4 and 4, the mean particle positions are stationary, with $\overline{x}\approx x$, and the coarsegrained computation of the intrinsic diffusion coefficient $cg\mathrm{\_}D$ is highly accurate, even if the coarsegrained concentration is not yet stationary. The mean diffusion coefficient $\overline{D}$, the mean velocity $\overline{u}$, and the corresponding standard deviations are calculated from sequences of $cg\mathrm{\_}D$ and $u$ values at the ${N}_{a}$ points inside the spatial averaging interval averaged over the three sampling times. The results presented in Table 1 are highly accurate for both $cg\mathrm{\_}D$ and $u$.
$\overline{D}$  $\overline{u}$  

$\mathrm{\Delta}x$  BGRW  GRW  BGRW  GRW 
5.00e03  0.0001$\pm $ 2.25e19  0.0001$\pm $ 2.51e19  1.00$\pm $ 0.00  1.00$\pm $ 0.00 
2.50e03  0.0001$\pm $ 5.07e19  0.0001$\pm $ 6.21e19  1.00$\pm $ 0.00  1.00$\pm $ 0.00 
1.25e03  0.0001$\pm $ 1.16e18  0.0001$\pm $ 1.30e18  1.00$\pm $ 2.36e16  1.00$\pm $ 0.00 
6.25e04  0.0001$\pm $ 1.42e18  0.0001$\pm $ 1.78e18  1.00$\pm $ 2.36e16  1.00$\pm $ 0.00 
Since the implementation of the CGST approach with the BGRW microscopic description is simpler and free of overshooting errors in simulations of biodegradation reactions [Suciu and Radu,, 2022, Sect. 5.2], it will be used in the following to illustrate the spacetime upscaling in reactive transport problems.
4 Spacetime upscaling in onedimensional reactive transport problems
We consider the onedimensional advectiondiffusion transport in saturated porous media of two reactive chemical species governed at the finegrained scale by the coupled system of equations
$${\partial}_{t}{c}_{\nu}+{\partial}_{x}({c}_{\nu}u)D{\partial}_{x}^{2}{c}_{\nu}={R}_{\nu}({c}_{1},{c}_{2}),\nu =1,2.$$  (13) 
As shown in [Destouni and Graham,, 1997], the spacetime scales of the experimental measurements are determined by the observation method. The application of the CGST approach in modeling real experiments is the subject of forthcoming studies. In the following, we illustrate the CGST spacetime upscaling by hypothetical experiments with independent spatial and temporal scales which ensure a good resolution of the averages and are bellow the spacetime scales on which the reacting species concentrations undergo significant changes. Precisely, in the examples presented below, we choose the spatial scale $a$ as 1/30 of the length of the domain in the mean flow direction and the temporal scale $\tau $ as 1/11 of the total observation time.
The discrepancy between volume and CGST averages, ${\overline{c}}_{\nu}$ and $\u27e8{1}_{\nu}\u27e9$, at a given time point is quantified by relative errors computed with vector ${l}^{2}$norms as ${e}_{{c}_{\nu}}=\Vert {\overline{c}}_{\nu}\u27e8{1}_{\nu}\u27e9\Vert /\Vert \u27e8{1}_{\nu}\u27e9\Vert $, as well as by maximum relative differences
$${\epsilon}_{{c}_{\nu}}=\mathrm{max}({\overline{c}}_{\nu}\u27e8{1}_{\nu}\u27e9)/\u27e8{1}_{\nu}\u27e9(\mathrm{arg}(\mathrm{max}(\u27e8{1}_{\nu}\u27e9{\overline{c}}_{\nu}))),\nu =1,2.$$ 
The influence of $a$ and $\tau $ on discrepancy is analyzed for a specific problem in B.
4.1 Bimolecular reaction with conservation of the total mass
Let us begin by considering a simple reaction where the first species is consumed and the second species is augmented at the same rate. The corresponding reaction rates are given by ${R}_{1}={K}_{r}{c}_{1}{c}_{2}^{{\ast}^{2}}$ and ${R}_{2}={K}_{r}{c}_{1}{c}_{2}^{{\ast}^{2}}$, with constant ${K}_{r}$. Summing up the two equations (13) one obtains an advectiondiffusion equation for the passive transport of the sum $({c}_{1}+{c}_{2})$, that is, the total mass is conserved.
Similarly to the derivation of Eq. (3), we find that the CGST averages of the two molecular species satisfy the particular form of Eq. (2),
$${\partial}_{t}\u27e8{1}_{\nu}\u27e9+{\partial}_{x}(\u27e8{1}_{\nu}\u27e9\overline{\xi})=\delta {1}_{\nu},\nu =1,2,$$  (14) 
where we use the notation ${1}_{\nu}$ to specify the molecular species when integrating particle trajectories in Eq. (5) to compute the CGST average concentrations. According to (1), the terms $\delta {1}_{\nu}$ account for creation and consumption of particles, hence, they can be interpreted as spacetime upscaled reaction terms. We remark that while the finegrained concentrations ${c}_{\nu}$ are governed by the advectiondiffusionreaction equations (13), the balance equations (14) verified by the CGST concentrations $\u27e8{1}_{\nu}\u27e9$ have the general form of an advectionreaction equation.
For consistency reasons, the conservation of the sum $({c}_{1}+{c}_{2})$ of the finegrained concentrations should imply the conservation of the upscaled concentrations. The latter implies $\delta {1}_{1}+\delta {1}_{2}=0$ and the passive transport of the sum $(\u27e8{1}_{1}\u27e9+\u27e8{1}_{2}\u27e9)$. To verify this statement, we compute CGST averages on the BGRW description associated to Eqs. (13), with constant parameters $u=1$, $D=0.01$, and ${K}_{r}=0.1$. The two molecular species are represented, as in Section 3.3 above, by one mole of particles uniformly distributed on the onedimensional lattice at the initial time. As boundary conditions we fix the concentration on $l=1$ to its initial value and set the no flux condition ${n}_{\mathcal{L},k}={n}_{\mathcal{L}1,k}$ on the outflow boundary on $l=\mathcal{L}$. We use the domain $\mathrm{\Omega}=(0,1)$, total time $T=1.1$, spatial scale $a=0.03$, and temporal scale $\tau =0.1$. The reactive transport is solved by a splitting scheme which alternates BGRW transport steps with reaction steps computed deterministically (a onedimensional version of the twodimensional splitting approach proposed in [Suciu and Radu,, 2022, Sect., 3.1]). The computations are carried out with $\mathrm{\Delta}x=2.50\cdot {10}^{3}$ and $\mathrm{\Delta}t=1.56\cdot {10}^{4}$.
Since the concentration is not vanishing in the domain during the computation, the intrinsic diffusion coefficient and the velocity can be calculated as in Section 3 above. One obtains accurate estimates $D=9.96\cdot {10}^{3}\pm 1.42\cdot {10}^{3}$ and $u=9.98\cdot {10}^{1}\pm 1.70\cdot {10}^{2}$, which certifies the correctness of the CGST averaging procedure. In Figs. 6 and 6, the normalized CGST averages $\u27e8{1}_{\nu}\u27e9/\mathcal{N}$ are compared to the spatial averages computed, similarly to (11), by ${\overline{c}}_{\nu}=(1/2a){\sum}_{i=1}^{{N}_{a}}{n}_{i}$, at the center $t$ of the temporal interval. One remarks that the conservation of the sum of the two upscaled species concentrations is verified with high accuracy. As shown in Table 2, the volume and the CGST averages are practically indistinguishable, with relative differences of the order of ${10}^{9}$. Hence, for small temporal variations of the concentration, as in the example presented here, the volume average alone is an appropriate model of the experimental measurements.
$t$  ${e}_{{c}_{1}}$  ${\epsilon}_{{c}_{1}}$  ${e}_{{c}_{2}}$  ${\epsilon}_{{c}_{2}}$ 

0.1  2.44e09  7.54e09  2.44e09  7.54e09 
0.5  1.81e09  3.99e09  1.81e09  3.99e09 
1  7.54e10  2.53e09  7.54e10  2.53e09 
4.2 Reactive transport with Monod reactions
In the following, we apply the CGST upscaling procedure to aerobic biodegradation processes in the subsurface [Wiedemeier et al.,, 1999]. An electron donor contaminant (e.g., benzene) of concentration ${c}_{1}$ is consumed by the biomass contained in the aqueous system, hereafter assumed constant in time and uniformly distributed in $\mathrm{\Omega}$ with ${c}_{bio}=1$, in the presence of an electron acceptor (oxygen) of concentration ${c}_{2}$. The electron donor and acceptor are consumed at rates ${R}_{1}=\theta {\alpha}_{1}\mu $ and ${R}_{2}=\theta {\alpha}_{2}\mu $, where $\theta $ is the volumetric water content, ${\alpha}_{1}$ and ${\alpha}_{2}$ are stoichiometric constants, and
$$\mu =\frac{{c}_{1}}{{M}_{1}+{c}_{1}}\frac{{c}_{2}}{{M}_{2}+{c}_{2}}$$  (15) 
is a Monod term with parameters ${M}_{1}$ and ${M}_{2}$. For the purpose of the present illustration of the CGST procedure, we consider the example used in [Bause and Knabner,, 2004, Brunner et al.,, 2012] and we set the parameters of the reaction system to ${\alpha}_{1}=5$, ${\alpha}_{2}=0.5$, ${M}_{1}={M}_{2}=0.1$, constant $\theta =1$ for saturated flow, and variable $\theta $ for partially saturated soils.
4.2.1 Monod reactions in saturated flow regime
The onedimensional biodegradation problem presented below is formulated with the same constant transport parameters $u$ and $D$, computational domain $\mathrm{\Omega}$, final time $T$, spatial and temporal scales $a$ and $\tau $, and is solved with the same numerical scheme as in the previous section. The initial contaminant concentration is set to one mole uniformly distributed on the lattice sites $l=1,\mathrm{\dots},\mathrm{\Delta}\mathcal{L}$, with $\mathrm{\Delta}\mathcal{L}=(\mathcal{L}1)/10$, $\mathcal{L}=1/\mathrm{\Delta}x+1$, and zero otherwise. Complementary, the oxygen concentration is set to one mole uniformly distributed on the lattice sites $l=\mathrm{\Delta}\mathcal{L}+1,\mathrm{\dots},\mathcal{L}$, and zero otherwise. The condition ${n}_{\mathcal{L},k}={n}_{\mathcal{L}1,k}$ is imposed on the outflow boundary and the first $\mathrm{\Delta}\mathcal{L}$ sites are set to the initial condition after each time step $k$. One simulates in this way a constant injection of contaminant near the inflow boundary, the consumption of the contaminant and oxygen through the biodegradation process in $\mathrm{\Omega}$, and the evacuation of the mixture of contaminant and oxygen at the outflow boundary. This is a simplified onedimensional version of the twodimensional biodegradation problem for saturated aquifers presented in [Suciu and Radu,, 2022, Sect. 5.2].
The coarsegrained averages are compared to the volume averages in Figs. 8 and 8. Differences between the two averages are visible mainly at the first sampling time. The two averages are close to each other in regions without reactions between contaminant and oxygen, e.g., contaminant concentration for $$ at $t=1$ (Fig. 8) and oxygen concentration for $x>0.4$ at $t=0.1$ (Fig. 8). A more precise evaluation of the discrepancy is obtained by computing relative differences. Table 3 shows that while the global differences ${e}_{{c}_{\nu}}$, computed with ${l}^{2}$norms are, with a single exception, smaller than 10%, the maximum relative differences ${\epsilon}_{{c}_{\nu}}$ are overall larger, with values between 17% and 57%. Similar results are obtained in one and twodimensional simulations of the biodegradation in heterogeneous aquifers presented in the extended version of this study [Suciu et al., 2023b, ]. Thus, even though the flow velocity and the dispersion coefficient are constant, timedependent reactions make the transport process nonstationary and the volume average alone could be an inaccurate model of the measurement.
$t$  ${e}_{{c}_{1}}$  ${\epsilon}_{{c}_{1}}$  ${e}_{{c}_{2}}$  ${\epsilon}_{{c}_{2}}$ 

2  0.0992  0.1753  0.0634  0.5716 
4  0.0367  0.4012  0.0400  0.2435 
6  0.0199  0.3572  0.3330  0.2574 
4.2.2 Monod reactions in variably saturated soils
One expects that for timedependent reactions in non stationary flow regimes the discrepancy between volume and CGST averages will be even more pronounced. As an example, we consider a onedimensional model problem for biodegradation processes in soil columns. The water flow in variably saturated porous media is described by the onedimensional Richards equation
$${\partial}_{t}\theta (\psi ){\partial}_{z}\left[K(\theta (\psi )){\partial}_{z}(\psi +z)\right]=0,$$  (16) 
where $\psi (z,t)$ is the pressure head expressed in length units, $\theta $ is the volumetric water content, $K$ stands for the hydraulic conductivity of the medium, and $z$ is the height oriented positively upward. The water flux given by Darcy’s law is $q=K(\theta (\psi ))\frac{\partial}{\partial z}(\psi +z)$. The reactive transport is now governed by Richards equation (16) coupled with the transport equations
$${\partial}_{t}\left[\theta (\psi ){c}_{\nu}\right]+{\partial}_{z}(q{c}_{\nu})D{\partial}_{z}^{2}{c}_{\nu}={R}_{\nu}({c}_{1},{c}_{2},\theta ),\nu =1,2.$$  (17) 
We consider in the following the relationships defining the water content $\theta (\psi )$ and the hydraulic conductivity $K(\theta (\psi ))$ given by the van GenuchtenMualem model
$$  (18) 
$$  (19) 
where $\mathrm{\Theta}=(\theta {\theta}_{res})/({\theta}_{sat}{\theta}_{res})$ is the normalized water content, ${\theta}_{res}$ and ${\theta}_{sat}$ are the residual and the saturated water content, ${K}_{sat}$ is the hydraulic conductivity of the saturated soil, and $\alpha $, $n$, $m=11/n$ are model parameters depending on the soil type. In this example, we consider the parameters of a silt loam used in previously published studies [List and Radu,, 2016, Illiano et al.,, 2020, Suciu et al.,, 2021, Suciu and Radu,, 2022], ${\theta}_{sat}=0.396$, ${\theta}_{res}=0.131$, $\alpha =0.423$, $n=2.06$, ${K}_{sat}=4.96\cdot {10}^{2}$. The soil heterogeneity is modeled by a fixed realization of a random saturated hydraulic conductivity with mean equal to the parameter ${K}_{sat}$ of the loam soil model. The realization is computed with the Kraichnan method described in [Suciu,, 2019, Appendix C.3.1.2], for a random $\mathrm{ln}K$ field with Gaussian correlation specified by the correlation length $\lambda =0.1$ and the variance ${\sigma}^{2}=0.5$. The dispersion is parameterized by the constant coefficient $D=0.001$. The system of coupled equations (1617), parameterized by the model (1819), with reactions governed by the Monod model (15) using the same parameters as in the saturated case presented above, is solved in the domain $\mathrm{\Omega}=(0,3)$ for the final time $T=6.6$. The spacetime scale of the CGST procedure is specified by $a=0.1$ and $\tau =0.6$.
The infiltration of the contaminated water in the column is driven by the variable pressure on the top boundary given by $\psi (3,t\le {t}_{1})=3+3.2t/{t}_{1}$ and $\psi (3,t>{t}_{1})=0.2$, with ${t}_{1}=T/3$. The initial pressure corresponds to the unsaturated profile $\psi (z,0)=z$ and the constant pressure condition $\psi (0,t)=\psi (0,0)$ is set on the bottom boundary. The initial contaminant concentration is set to one mole uniformly distributed on the lattice sites $l=\mathcal{L}\mathrm{\Delta}\mathcal{L},\mathrm{\dots},\mathcal{L}$, with $\mathrm{\Delta}\mathcal{L}=(\mathcal{L}1)/10$, $\mathcal{L}=3/\mathrm{\Delta}z+1$, on the top of the domain $\mathrm{\Omega}$ and zero otherwise. The oxygen concentration is set to one mole uniformly distributed on the lattice sites $l=1,\mathrm{\dots},\mathcal{L}\mathrm{\Delta}\mathcal{L}1$ and zero otherwise. The top $\mathrm{\Delta}\mathcal{L}$ sites are set to the initial condition after each time step $k$ and the condition ${n}_{1,k}={n}_{2,k}$ is imposed on the bottom boundary for both species. With these, one simulates a continuous injection of contaminant over a length $\mathrm{\Delta}\mathcal{L}\mathrm{\Delta}z$ on the top of the column, assuming that the incoming contaminated water does not contain oxygen, in conditions of free drainage at the bottom of the column.
The flow equation (16) is solved with the iterative GRW algorithm for onedimensional flows, implemented as an $L$scheme linearization approach (see [Pop et al.,, 2004, List and Radu,, 2016]) capable to handle the nonlinearity and the degeneracy of the Richards equation [Suciu et al.,, 2021, Sect. 2]. At each time step $k$, the solution $(\psi ,\theta )$ of the Richards equation is transmitted to the reactive transport solver. The latter is an iterative $L$scheme, which at each iteration step solves the system (17) with an operator splitting procedure consisting of a BGRW solution of the advectiondiffusion process and a deterministic solution of the reaction step (this is the onedimensional version of the twodimensional iterative scheme introduced in [Suciu and Radu,, 2022, Sect. 3.2]). The spatial step is fixed to $\mathrm{\Delta}z=1.25\cdot {10}^{2}$. The time steps, determined as the minimum of the time steps required by the GRW flow solver and the BGRW transport solver lay between $\mathrm{\Delta}t=7.42\cdot {10}^{3}$ and $\mathrm{\Delta}t=7.3\cdot {10}^{3}$. The number of iterations required for the convergence of the two $L$schemes ranges between 100 and 450.
The flow solutions and the spacetime upscaled concentrations are presented in Figs. 1012. The pressure profiles recorded at the three sampling times shown in Fig. 10 capture the spatial heterogeneity of the flow solution and the transition from the unsaturated regime ($$) to the saturated regime ($\psi \ge 0$). The resulting flow velocity (Fig. 10) is strongly variable in both space and time. The complexity of the unsaturated/saturated flow impacts on the distribution of the reacting species (see Figs. 12 and 12), which inherits the heterogeneity and the spatiotemporal variability of the flow. In these conditions, the discrepancy between the volume and the CGST average is dramatically large at all the tree sampling times, with maximum relative differences between 59% and 99% (Table 4). Thus, averaging over the spatial interval $(xa,x+a)$ at given time $t$ is completely irrelevant for the observation/measurement made at the spatiotemporal scale $(a,\tau )$ centered at $(x,t)$. Excepting special situations, such as the case of almost homogeneous reactions with slow time variations presented in Figs. 6 and 6, the time average has to be considered as well in onedimensional modeling of experimental observations.
$t$  ${e}_{{c}_{1}}$  ${\epsilon}_{{c}_{1}}$  ${e}_{{c}_{2}}$  ${\epsilon}_{{c}_{2}}$ 

2  0.5054  0.9895  0.1127  0.5867 
4  0.3420  0.6950  0.2392  0.6057 
6  0.3136  0.7702  0.3745  0.6576 
5 CGST averaging in twodimensional transport problems
The twodimensional model for reactive transport with biodegradation reactions, similar to the onedimensional model considered above, is governed by the equations
${\partial}_{t}\theta (\psi )\nabla \cdot [K(\theta (\psi )\nabla (\psi +z)]=0,$  (20)  
${\partial}_{t}\left[\theta (\psi ){c}_{\nu}\right]\nabla \cdot (D\nabla {c}_{\nu}\mathbf{q}{c}_{\nu})={R}_{\nu}({c}_{1},{c}_{2},\theta ),\mathbf{q}=K(\theta (\psi )\nabla (\psi +z),\nu =1,2.$  (21) 
Equations (2021) are parameterized with the van GenuchtenMualem model (1819) and the same Monod reaction system as in the previous sections. The twodimensional CGST averaging procedure, first verified by computing the velocity components and the diffusion coefficients, will be used in the following to compute spacetime upscaled concentrations.
5.1 Verification of the two dimensional CGST averaging procedure
The twodimensional CGST averaging procedure corresponds to $d=2$ in the general definition (1). Similarly to the derivation of the onedimensional stationary diffusion equation (8) given in Section 3, one obtains the diffusion equation
$${\partial}_{{x}_{\alpha}}{\partial}_{{x}_{\beta}}({D}_{\alpha \beta}c)=0,$$ 
with diffusion tensor defined by ${D}_{\alpha \beta}=\u27e8{x}_{\alpha}{\xi}_{\beta}\u27e9/\u27e81\u27e9$. The numerical implementation of the CGST averaging procedure is a straightforward extension of that used in the onedimensional case.
We perform BGRW/GRW simulations in the twodimensional domain $\mathrm{\Omega}=(0,1)\times (0,1)$ of the advectiondiffusion process defined by constant diagonal diffusion tensor with components ${D}_{11}={D}_{22}=0.0001$ and constant velocity $(u,v)=(1,0)$. The BGRW/GRW simulations are done with the same number of particles and spacetime discretization and the CGST averages are computed with the same parameters $a=0.05$ and $\tau =0.1$ as in the one dimensional case. The components of the diffusion tensor (Tables 5 and 6) and of the velocity (Table 7) computed with the CGST procedure are, again, in very good agreement with the nominal coefficients of the advectiondiffusion process.
$\mathrm{\Delta}x$  $\overline{{D}_{11}}$  $\overline{{D}_{12}}$  $\overline{{D}_{21}}$  $\overline{{D}_{22}}$ 

5.00e03  0.0001$\pm $ 1.59e19  1.59e20$\pm $ 8.13e20  0.00$\pm $ 0.00  0.0001$\pm $ 0.00 
2.50e03  0.0001$\pm $ 1.31e19  1.06e19$\pm $ 2.92e19  0.00$\pm $ 0.00  0.0001$\pm $ 0.00 
1.25e03  0.0001$\pm $ 3.02e19  5.93e20$\pm $ 2.36e19  0.00$\pm $ 0.00  0.0001$\pm $ 2.87e20 
6.25e04  0.0001$\pm $ 7.47e19  1.23e19$\pm $ 3.02e19  0.00$\pm $ 0.00  0.0001$\pm $ 1.44e20 
$\mathrm{\Delta}x$  $\overline{{D}_{11}}$  $\overline{{D}_{12}}$  $\overline{{D}_{21}}$  $\overline{{D}_{22}}$ 

5.00e03  0.0001$\pm $ 1.20e19  3.57e21$\pm $ 7.29e20  0.00$\pm $ 0.00  0.0001$\pm $ 1.44e20 
2.50e03  0.0001$\pm $ 1.63e19  5.21e20$\pm $ 2.63e19  0.00$\pm $ 0.00  0.0001$\pm $ 2.87e20 
1.25e03  0.0001$\pm $ 2.46e19  1.09e19$\pm $ 2.73e19  0.00$\pm $ 0.00  0.0001$\pm $ 1.44e20 
6.25e04  0.0001$\pm $ 9.66e19  2.07e19$\pm $ 3.56e19  0.00$\pm $ 0.00  0.0001$\pm $ 1.44e20 
$\overline{u}$  $\overline{v}$  

$\mathrm{\Delta}x$  BGRW  GRW  BGRW  GRW 
5.00e03  1.00$\pm $ 0.00  1.00$\pm $ 0.00  2.71e19$\pm $ 0.00  0.00$\pm $ 0.00 
2.50e03  1.00$\pm $ 1.18e16  1.00$\pm $ 0.00  1.39e19$\pm $ 0.00  0.00$\pm $ 0.00 
1.25e03  1.00$\pm $ 0.00  1.00$\pm $ 0.00  1.14e19$\pm $ 2.55e35  0.00$\pm $ 0.00 
6.25e04  1.00$\pm $ 2.36e16  1.00$\pm $ 0.00  9.55e20$\pm $ 1.28e35  0.00$\pm $ 0.00 
5.2 Application to reactive transport in variably saturated soils
To illustrate the twodimensional CGST averaging procedure, we simulate the same biodegradation process as in the onedimensional case in the domain $\mathrm{\Omega}=(0,2)\times (0,3)$. We consider the same soil model, vertical profile of the pressure, the same top and bottom pressure boundary conditions, and noflow conditions on the left and right boundaries. The initial concentrations for contaminant and oxygen in the rectangle ${\mathrm{\Omega}}_{1}=[0.5,1.5]\times [2.7,3]$ and in the surrounding region take the same values as on the top $\mathrm{\Delta}\mathcal{L}$ sites and on the remaining sites of the lattice in the onedimensional problem and concentrations in ${\mathrm{\Omega}}_{1}$ are set to the initial condition after each time step. Noflux conditions are imposed on the bottom and left/right side boundaries for both the contaminant and the oxygen concentrations. CGST averages are computed with $a=0.1$ and $\tau =12$, for a total simulation time $T=132$. The flow and the BGRW solutions are computed with the flow $L$scheme from [Suciu et al.,, 2021, Sect. 4.1] coupled to the iterative splitting procedure introduced in [Suciu and Radu,, 2022, Sect. 3.2].
CGST averages are computed along a vertical sampling line placed at the center of the horizontal edge, as well as along a vertical line close to the right boundary of the domain. Figures 13  15 show the flow solution and the finegrained contaminant and oxygen concentrations at the final time $T=132$. Significant differences between the CGST and volume averages at the three sampling times can be observed for both contaminant and oxygen concentrations presented in Figs. 1719 in regions where the two species interact and negligible differences outside the reaction zone (indicated by constant concentrations). One remarks also that in case of decentered sampling, near the right boundary (Figs. 19 and 19), these differences are generally larger than for centered sampling line (Figs. 17 and 17). Table 8 quantifies these visual observations and warns that disregarding the time average when modeling measurements through twodimensional simulations of reactive transport in partially saturated soils may result in discrepancies up to 80% with respect to the CGST averaging approach.
centered  decentered  

$t$  ${e}_{{c}_{1}}$  ${\epsilon}_{{c}_{1}}$  ${e}_{{c}_{2}}$  ${\epsilon}_{{c}_{2}}$  ${e}_{{c}_{1}}$  ${\epsilon}_{{c}_{1}}$  ${e}_{{c}_{2}}$  ${\epsilon}_{{c}_{2}}$ 
40  0.3466  0.8041  0.0675  0.1389  0.6080  0.7670  0.1320  0.5471 
80  0.1198  0.1636  0.0536  0.3479  0.1426  0.1737  0.0776  0.1354 
120  0.0654  0.0880  0.0280  0.0462  0.0862  0.1062  0.0440  0.2204 
6 Summary
The CGSTGRW approach proposed in this article performs the spacetime upscaling adapted to the spatiotemporal scale of a hypothetical measurement. The CGST averages over symmetrical volumes and time intervals of the microscopic description of a corpuscular system given by piecewise analytic time functions are a.e. continuous fields which verify identities similar to the local balance equations of the fluid dynamics. If the averages correspond to upscaled concentrations, these identities allow us to derive theoretical expressions for the flow velocity and the intrinsic diffusion coefficients of the system of particles. The microscopic description of reactive transport in porous media consists of GRW simulations of the advectivediffusive movement of the actual number of molecules involved in reactions. The chemical reactions are modeled deterministically by solving the reaction system with concentrations expressed as mole fractions of the number of reacting molecules. The main results of this study are summarized as follows.

1.
The codes for the CGSTGRW numerical modeling are verified by computing the theoretical expressions for flow velocity and intrinsic diffusion coefficients and one finds an excellent agreement with the nominal values of these input parameters.

2.
If the reacting system consists of biomolecular reactions in a uniform flow, the initial distribution of the reacting species is uniform, and the boundary conditions force constant concentration solutions, the CGST average is practically indistinguishable from the classical volume average.

3.
For biodegradation reactions in saturated flows, the CGST averages may significantly differ from the volume averages even if the flow is stationary, due to the timedependence of the reaction process.

4.
The discrepancy between volume and CGST averages may be dramatically increased if the biodegradation process takes place in partially saturated soils with timedependent flow velocity.

5.
According to (4), the CGST average concentration is equivalent to a spacetime average over the continuous field of the finegrained concentration. This renders possible comparisons with time averages of classical volume sampling solutions obtained by integrating continuous fields or with time averages over solutions provided by homogenization methods. Kernel density estimates of concentrations in particle tracking simulations averaged over symmetrical time intervals can be compared as well with CGST averages computed with the same spatial and temporal scales.
Appendix A CGST averaging on $d$dimensional cubes
We consider the evolution of a system of $\mathcal{N}$ particles during the time interval $I=[0,T]$. Let ${I}_{i}=[{t}_{i}^{+},{t}_{i}^{}]$, $0\le {t}_{i}^{+}\le {t}_{i}^{}$, $i=1,\mathrm{\dots},\mathcal{N}$, be the interval of existence of each particle, where ${t}_{i}^{+}$ and ${t}_{i}^{}$ are the birth and dead times, respectively. These time points could be, for instance, those of generation and consumption of a molecule in chemical reactions. A kinematic microscopic description of a physical system can be achieved by arbitrary time functions ${\phi}_{i}(t):I\mapsto \mathbb{R}$, with ${\phi}_{i}(t)=0$ $\forall t\in I\setminus {I}_{i}$, $i=1,\mathrm{\dots},\mathcal{N}$. Particular cases of functions ${\phi}_{i}$ are the components of the particle’s trajectory, ${x}_{\alpha i}:I\mapsto \mathbb{R}$, and the corresponding velocity components, ${\xi}_{\alpha i}:I\mapsto \mathbb{R}$, ${\xi}_{\alpha i}(t)=d{x}_{\alpha i}(t)/dt$. The existence intervals ${I}_{i}$ may contain finite number of points consisting of sets of zero Lebesgue measure ${I}_{i}^{\ast}=\{{t}_{i1},\mathrm{\dots},{t}_{ik},{t}_{i(k+1)},\mathrm{\dots}{t}_{in}\}\subset {I}_{i}$ where the functions ${\phi}_{i}$ undergo discontinuous variations. For instance, this is the case of velocity changes at successive time steps in random walk simulations. The restrictions of ${\phi}_{i}$ to intervals $({t}_{ik},{t}_{i(k+1)})$ are assumed to be analytic functions. So, we consider a kinematic description given by sets of piecewise analytic functions associated to the system of particles.
Let $$ be an open $d$dimensional cube with side $2a$, centered at $\mathbf{x}$. Similarly to the approach for coarsegrained averaging over phase points describing ensembles in statistical mechanics, for the case ${I}_{i}=I,\forall i=1,\mathrm{\dots},\mathcal{N}$, presented in [Vamoş,, 2007, Chap. 5], we define the CGST average over the cube $C(\mathbf{x},a)$ and the time interval $(t\tau ,t+\tau )$, $$,
$$\u27e8\phi \u27e9(\mathbf{x},t;a,\tau )=\frac{1}{2\tau {(2a)}^{d}}\sum _{i=1}^{\mathcal{N}}{\int}_{t\tau}^{t+\tau}{\phi}_{i}({t}^{\prime})\prod _{\alpha =1}^{d}{\chi}_{\alpha}({x}_{\alpha i}({t}^{\prime}))d{t}^{\prime},$$  (22) 
where ${\chi}_{\alpha}$ is the characteristic function of the interval $({x}_{\alpha}a,{x}_{\alpha}+a)$.
A particle enters (leaves) the cube $C(\mathbf{x},a)$ at time instants ${u}_{i}$ which solve the system of equations
$${x}_{\alpha i}({u}_{i}){x}_{\alpha}\pm a=0,\alpha =1,2,\mathrm{\dots},d.$$  (23) 
We denote by ${U}_{i}^{\prime}=\{{u}_{i1}^{\prime},{u}_{i2}^{\prime},\mathrm{\dots},{u}_{i{n}^{\prime}}^{\prime}\}$ and ${U}_{i}^{\prime \prime}=\{{u}_{i1}^{\prime \prime},{u}_{i2}^{\prime \prime},\mathrm{\dots},{u}_{i{n}^{\prime \prime}}^{\prime \prime}\}$ the sets of time points corresponding to instances at which the $i$th particles enters or leaves the cube $C(\mathbf{x},a)$. It follows that, for a fixed $\mathbf{x}$, the integrand in Eq. (22),
$${G}_{i}(\mathbf{x},{t}^{\prime})={\phi}_{i}({t}^{\prime})\prod _{\alpha =1}^{d}{\chi}_{\alpha}({x}_{\alpha i}({t}^{\prime}))={\phi}_{i}({t}^{\prime}){\chi}_{C(\mathbf{x},a)}({\mathbf{x}}_{i}({t}^{\prime})),$$  (24) 
is a continuous time function in $(t\tau ,t+\tau )$, except at a finite number of points $\{{t}_{i}^{+},{t}_{i}^{}\}\bigcup {U}_{i}^{\prime}\bigcup {U}_{i}^{\prime \prime}\bigcup {I}_{i}^{\ast}$ where it has finite jump discontinuities. Hence, $G(\mathbf{x},{t}^{\prime})$ is a function with bounded variation and Riemann integrable.
Using Eq. (22) and Eq. (24), the time derivative of the CGST average can be expressed as
$${\partial}_{t}\u27e8\phi \u27e9(\mathbf{x},t;a,\tau )=\frac{1}{2\tau {(2a)}^{d}}\sum _{i=1}^{\mathcal{N}}\left[{G}_{i}(\mathbf{x},t+\tau ){G}_{i}(\mathbf{x},t\tau )\right].$$  (25) 
The continuity of the time derivative (25) is determined by that of the function ${G}_{i}$. According to (24), ${G}_{i}$ is not continuous when ${\phi}_{i}$ is discontinuous and ${\chi}_{C(\mathbf{x},a)}$ is nonvanishing or, conversely, ${\chi}_{C(\mathbf{x},a)}$ is discontinuous and ${\phi}_{i}$ is nonvanishing. The first situation occurs if the particle is created or consumed at $t={t}_{i}^{\pm}$ or undergoes discontinuous variations at times $t\in {I}_{i}^{\ast}$ inside the cube $C(\mathbf{x},a)$. The second situation occurs during the time intervals $\subset {I}_{i}\setminus {I}_{i}^{\ast}$ when the particle lies on the surface of the cube, $\partial C(\mathbf{x},a)$. Taking into account the time arguments of ${G}_{i}$ in (25), ${\partial}_{t}\u27e8\phi \u27e9$ is discontinuous on the set ${\mathrm{\Omega}}^{\prime}={\cup}_{i=1}^{\mathcal{N}}{\mathrm{\Omega}}_{i}^{\prime}$, where
${\mathrm{\Omega}}_{i}^{\prime}=$  $\{\{{t}_{i}^{\pm}\tau ,{t}_{i}^{\pm}+\tau \}\cap (\tau ,T\tau )\}\times C({\mathbf{x}}_{i}({t}_{i}^{\pm}),a)\}\cup $  
$\{\{\{{t}_{ik}\tau ,{t}_{ik}+\tau \}{t}_{ik}\in {I}_{i}^{\ast}\}\cap (\tau ,T\tau )\}\times C({\mathbf{x}}_{i}({t}_{i}^{\ast}),a)\}\cup $  
$\{\{[{t}_{i}^{+}\pm \tau ,{t}_{i}^{}\pm \tau ]\cap (\tau ,T\tau )\}\times \partial C({\mathbf{x}}_{i}(t\mp \tau ),a)\}\cup $  
$\{\{\{[{t}_{i}^{+}\pm \tau ,{t}_{i1}\pm \tau )\cup ({t}_{i1}\pm \tau ,{t}_{i2}\pm \tau )\cup \mathrm{\dots}$  
$\cup ({t}_{in}\pm \tau ,{t}_{i}^{}\pm \tau ]\}\cap (\tau ,T\tau )\}\times \partial C({\mathbf{x}}_{i}(t\mp \tau ),a)\}.$ 
Since ${\mathrm{\Omega}}^{\prime}$ has zero Lebesgue measure in ${\mathbb{R}}^{d}\times (\tau ,T\tau )$, the time derivative ${\partial}_{t}\u27e8\phi \u27e9$ is a.e. continuous.
As a function with bounded variation, ${G}_{i}$ may be decomposed as a sum of a continuous function and a jump function, ${G}_{i}={G}_{i}^{\prime}+{G}_{i}^{\prime \prime}$. Accordingly,
$${\partial}_{t}\u27e8\phi \u27e9={({\partial}_{t}\u27e8\phi \u27e9)}^{\prime}+{({\partial}_{t}\u27e8\phi \u27e9)}^{\prime \prime}.$$  (26) 
Since, except at discontinuity points, ${G}_{i}$ is an analytic function, according to the fundamental theorem of Lebesgue integral calculus, the continuous part ${G}_{i}^{\prime}$ is absolutely continuous. Hence,
$${({\partial}_{t}\u27e8\phi \u27e9)}^{\prime}=\u27e8\frac{d\phi}{dt}\u27e9.$$  (27) 
The second term of Eq. (26) consists of positive contributions when particles enter into the cube $C(\mathbf{x},a)$ and negative contributions when they leave the cube, as well as contributions of other three types of discontinuities, which will be described below. According to Eq. (22), ${({\partial}_{t}\u27e8\phi \u27e9)}^{\prime \prime}$ can be written as
$${({\partial}_{t}\u27e8\phi \u27e9)}^{\prime \prime}=\frac{1}{2\tau {(2a)}^{d}}\sum _{i=1}^{\mathcal{N}}\sum _{\alpha =1}^{d}\left[\sum _{{u}_{\alpha}\in {W}_{i}^{\prime}}{\phi}_{i}({u}_{\alpha})\sum _{{u}_{\alpha}\in {W}_{i}^{\prime \prime}}{\phi}_{i}({u}_{\alpha})\right]+{({\partial}_{t}\u27e8\phi \u27e9)}_{gen}+{({\partial}_{t}\u27e8\phi \u27e9)}_{disc}+\u03f5,$$  (28) 
where ${W}_{i}^{\prime}={U}_{i}^{\prime}\bigcap (t\tau ,t+\tau )$ and ${W}_{i}^{\prime \prime}={U}_{i}^{\prime \prime}\bigcap (t\tau ,t+\tau )$. The term ${({\partial}_{t}\u27e8\phi \u27e9)}_{gen}$ accounts for generation/consumption of molecules in chemical reactions. It stems from variations of ${G}_{i}$ when the molecule is generated/consumed in the interior of $C(\mathbf{x},a)$ during the interval $(t\tau ,t+\tau )$:
$${\mathrm{\Delta}}^{+}{G}_{i}={\phi}_{i}({t}^{+}){\chi}_{{}_{C(\mathbf{x},a)}}({\mathbf{x}}_{i}({t}^{+})){\chi}_{{}_{(t\tau ,t+\tau )}}({t}^{+}),{\mathrm{\Delta}}^{}{G}_{i}={\phi}_{i}({t}^{}){\chi}_{{}_{C(\mathbf{x},a)}}({\mathbf{x}}_{i}({t}^{})){\chi}_{{}_{(t\tau ,t+\tau )}}({t}^{}),$$ 
where ${\chi}_{C(\mathbf{x},a)}$ and ${\chi}_{(t\tau ,t+\tau )}$ are the characteristic functions of the cube and of the timeaveraging interval, respectively. Hence,
$${({\partial}_{t}\u27e8\phi \u27e9)}_{gen}=\frac{1}{2\tau {(2a)}^{d}}\sum _{i=1}^{\mathcal{N}}\left[{\mathrm{\Delta}}^{+}{G}_{i}+{\mathrm{\Delta}}^{}{G}_{i}\right].$$ 
The term ${({\partial}_{t}\u27e8\phi \u27e9)}_{disc}$ describes variations due to the discontinuities of ${\phi}_{i}$ at points from ${I}_{i}^{\ast}$ [Vamoş et al.,, 1997, 2000]. This is, for instance, the case if ${\phi}_{i}={\xi}_{\alpha i}$ in random walk simulations when the velocity undergoes finite jumps at the end of each time step. In general,
$${({\partial}_{t}\u27e8\phi \u27e9)}_{disc}=\frac{1}{2\tau {(2a)}^{d}}\sum _{i=1}^{\mathcal{N}}\sum _{s\in {I}_{i}^{\ast}}\left[{\phi}_{i}(s0){\phi}_{i}(s+0)\right],$$ 
where ${\phi}_{i}(s0)$ and ${\phi}_{i}(s+0)$ are the limits from the left and from the right respectively. The term $\u03f5$ cancels spurious contributions ${\phi}_{i}({u}_{\alpha})$ which occur when the trajectory intersects an vertex or an edge of the $d$dimensional cube with $d\ge 2$. Since $\u03f5$ is defined on a set with zero Lebesgue measure, we assume that it is negligible. For instance in twodimensional random walk simulations, unless in the worst case when most of the particles are located in the four corners of the square, the relative number of spurious contributions, of the order $1/\mathcal{L}$, where $\mathcal{L}$ is the number of lattice points per side, is quite small and $\u03f5$ can be disregarded.
The CGST average defined by (22) depends on ${x}_{\alpha}$ through the time points ${u}_{i}^{\prime}$ and ${u}_{i}^{\prime \prime}$. The derivative of (23) w.r.t. ${x}_{\alpha}$ leads to $\partial {u}_{i}/\partial {x}_{\alpha}=1/(d{x}_{\alpha i}/d{u}_{i})=1/{\xi}_{\alpha i}$. If ${u}_{i}\in (t\tau ,t+\tau )$, then ${u}_{i}^{\prime}$ is the lower integration limit of the time integral in (22) and ${u}_{i}^{\prime \prime}$ is the upper limit. This implies the following expression of the partial space derivative,
$${\partial}_{{x}_{\alpha}}\u27e8\phi \u27e9=\frac{1}{2\tau {(2a)}^{d}}\sum _{i=1}^{\mathcal{N}}\left[\sum _{{u}_{\alpha}\in {W}_{i}^{\prime}}\frac{{\phi}_{i}({u}_{\alpha})}{{\xi}_{\alpha i}}+\sum _{{u}_{\alpha}\in {W}_{i}^{{}^{\prime \prime}}}\frac{{\phi}_{i}({u}_{\alpha})}{{\xi}_{\alpha i}}\right].$$  (29) 
As long as the particle’s velocity ${\xi}_{\alpha i}\ne 0$ is defined and the implicit function theorem can be applied to define $\partial {u}_{i}/\partial {x}_{\alpha}=1/{\xi}_{\alpha i}$, the spatial derivative ${\partial}_{{x}_{\alpha}}\u27e8\phi \u27e9$ exists and is continuous. These conditions are not met on sets ${\mathrm{\Omega}}_{i}^{\prime \prime}$ consisting of time intervals determined by $t\pm \tau $ with $t\in \{{t}_{i}^{+},{t}_{i}^{}\}\cup {I}^{\ast}$ and sets $\{\mathbf{x}\in \partial C({\mathbf{x}}_{i}(t),a)\}$ as well as in cases where particles move on the surface of the cube and Eq. (23) is identically satisfied. A detailed analysis and the construction of the sets ${\mathrm{\Omega}}_{i}^{\prime \prime}$ is rather involved [Vamoş,, 2007, Chap. 4] but, in the present context, it is not necessary. Indeed, since the surface $\partial C({\mathbf{x}}_{i}(t),a)$ has a null Lebesgue measure in ${\mathbb{R}}^{d}$, the union ${\mathrm{\Omega}}^{\prime \prime}={\cup}_{i=1}^{\mathcal{N}}{\mathrm{\Omega}}_{i}^{\prime \prime}$ also has a null Lebesgue measure in ${\mathbb{R}}^{d}\times (\tau ,T\tau )$ and the spatial derivatives ${\partial}_{{x}_{\alpha}}\u27e8\phi \u27e9$ are a.e. continuous. The continuity of the partial derivatives ${\partial}_{{x}_{\alpha}}\u27e8\phi \u27e9$ and ${\partial}_{t}\u27e8\phi \u27e9$ ensures the a.e. continuity of the CGST averages $\u27e8\phi \u27e9$.
From (29) and (28) one obtains
$${({\partial}_{t}\u27e8\phi \u27e9)}^{{}^{\prime \prime}}=\sum _{\alpha =1}^{d}{\partial}_{{x}_{\alpha}}\u27e8\phi {\xi}_{\alpha}\u27e9+{({\partial}_{t}\u27e8\phi \u27e9)}_{gen}+{({\partial}_{t}\u27e8\phi \u27e9)}_{disc}.$$  (30) 
By inserting (27) and (30) into (26), one obtains the following identity verified by the CGST averages,
$${\partial}_{t}\u27e8\phi \u27e9+\nabla \cdot \u27e8\phi \bm{\xi}\u27e9=\u27e8\frac{d\phi}{dt}\u27e9+\delta \phi ,$$  (31) 
where $\delta \phi ={({\partial}_{t}\u27e8\phi \u27e9)}_{gen}+{({\partial}_{t}\u27e8\phi \u27e9)}_{disc}$. Equation (31) is valid in ${\mathbb{R}}^{d}\times (\tau ,T\tau )\setminus ({\mathrm{\Omega}}^{\prime}\cup {\mathrm{\Omega}}^{\prime \prime})$, that is, a.e. in the spacetime domain on which the CGST averages $\u27e8\phi \u27e9$ are defined.
Appendix B Discrepancy between volume and CGST averages for different $a$ and $\tau $.
Discrepancies between the volume and CGST averages are computed for the onedimensional example of reactive transport governed by Monod model presented in Section 4.2.2. The influence of the temporal scale $\tau $ and of the spatial scale $a$ on the global relative difference ${e}_{{c}_{\nu}}$ is illustrated in the subplots from Fig. 1. Similarly, the influence of $\tau $ and $a$ on ${\epsilon}_{{c}_{\nu}}$ is illustrated in Fig. 2.
Figure Fig. 1 shows a monotone increase of ${e}_{{c}_{\nu}}$ with $\tau $ for fixed $a$ and a monotone decrease with $a$ for fixed $\tau $. One remarks larger ${e}_{{c}_{\nu}}$ values for contaminant ($\nu =1$) and smaller values for oxygen ($\nu =2$) at the first sampling time (that is, close to the contaminant source)
Figure Fig. 2 presents a rather irregular behavior of ${\epsilon}_{{c}_{\nu}}$ at the first sampling time. The maximum difference is less influenced by the increase of $\tau $ and $a$. This is especially notable for the contaminant, with ${\epsilon}_{{c}_{1}}$ values that remain almost constant with the decrease of $\tau $ and the increase of $a$.
Acknowledgements
The authors are grateful to three antonymous reviewers for their valuable comments. N. Suciu acknowledges the financial support of the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Grant SU 415/41 – 405338726 “Integrated global random walk model for reactive transport in groundwater adapted to measurement spatiotemporal scales”.
References
 Akagi and Oka, [2020] Akagi, G., Oka, T., 2020. Spacetime homogenization for nonlinear diffusion. Preprint arXiv:2007.09977. https://doi.org/10.48550/arXiv.2007.09977
 Aquino and Dentz, [2017] Aquino, T., Dentz, M., 2017. Chemical continuous time random walks. Physical review letters, 119(23), 230601. https://doi.org/10.1103/PhysRevLett.119.230601
 Auriault and Adler, [1995] Auriault, J.L., Adler, P.M., 1995. Taylor dispersion in porous media: analysis by multiple scale expansions. Advances in Water Resources, 18(4), 217226. https://doi.org/10.1016/03091708(95)000117
 Andričević, [1998] Andričević, R., 1998. Effects of local dispersion and sampling volume on the evolution of concentration fluctuations in aquifers. Water Resour. Res., 34(5), 1115–29. http://dx.doi.org/10.1029/98WR00260
 de Barros and Fiori, [2014] de Barros, F.P.J., Fiori, A., 2014. Firstorder based cumulative distribution function for solute concentration in heterogeneous aquifers: Theoretical analysis and implications for human health risk assessment. Water Resour. Res., 50(5), 4018–4037. http://dx.doi.org/10.1002/2013WR015024
 BayerRaich et al., [2004] BayerRaich, M., Jarsjö, J., Liedl, R., Ptak, T., Teutsch, G.2004. Average contaminant concentration and mass flow in aquifers from timedependent pumping well data: Analytical framework. Water Resour. Res., 40(8), W08303. http://dx.doi.org/10.1029/2004WR003095
 Bause and Knabner, [2004] Bause, M., Knabner, P., 2004. Numerical simulation of contaminant biodegradation by higher order methods and adaptive time stepping. Comput. Visual. Sci., 7(2), 61–78. https://doi.org/10.1007/s007910040139y
 Bensoussan et al., [2011] Bensoussan, A., Lions, J.L. and Papanicolaou, G., 2011. Asymptotic analysis for periodic structures (Vol. 374). American Mathematical Soc.
 Berkowitz, [2021] Berkowitz, B., 2022. HESS Opinions: Chemical transport modeling in subsurface hydrological systems – Space, time, and the “holy grail” of “upscaling”. Hydrol. Earth. Syst. Sci., 26, 21612180. https://doi.org/10.5194/hess2621612022
 Brunner et al., [2012] Brunner, F., Radu, F.A., Bause, M., Knabner, P., 2012. Optimal order convergence of a modified BDM1 mixed finite element scheme for reactive transport in porous media. Adv. Water Resour., 35, 163–171. http://dx.doi.org/10.1016/j.advwatres.2011.10.001
 Bringedal et al., [2016] Bringedal, C., Berre, I., Pop, I. S., Radu, F. A., 2016. Upscaling of nonisothermal reactive porous media flow under dominant Péclet number: the effect of changing porosity. Multiscale Model. Simul. 14(1), 502–533. https://doi.org/10.1137/15M1022781
 Caroni and Fiorotto, [2005] Caroni, E., Fiorotto, V., 2005. Analysis of concentration as sampled in natural aquifers. Transport Porous Media, 59(1),19–45. http://dx.doi.org/10.1007/s112420041119x
 Cushman, [1983] Cushman, J.H., 1983. Multiphase transport equations: I  general equation for macroscopic statistical, local spacetime homogeneity. Transport Theor. Stat. Phys., 12(1), 35–71. http://dx.doi.org/10.1080/00411458308212731
 Dabrowska and Rykala, [2021] Dabrowska, D., Rykala, W., 2021. A review of lysimeter experiments carried out on municipal landfill waste. Toxics, 9(2), 26. https://doi.org/10.3390/toxics9020026
 Destouni and Graham, [1997] Destouni, G., Graham, W., 1997. The influence of observation method on local concentration statistics in the subsurface. Water Resour .Res., 33(4),663–676. http://dx.doi.org/10.1029/96WR03955
 FernàndezGarcia and SanchezVila, [2011] FernàndezGarcia, D., SanchezVila, X., 2011. Optimal reconstruction of concentrations, gradients and reaction rates from particle distributions. J. Contam. Hydrol., 120, 99–114. https://doi.org/10.1016/j.jconhyd.2010.05.001
 Gray and Miller., [2013] Gray, W.G., Miller, C.T., 2013. A generalization of averaging theorems for porous medium analysis. Adv. Water Resour. 62, 227–37. https://doi.org/10.1016/j.advwatres.2013.06.006
 Gray et al., [2015] Gray, W.G., Dye, A.L., McClure, J.E., PyrakNolte, L.J., Miller, C.T., 2015. On the dynamics and kinematics of twofluidphase flow in porous media. Water Resour. Res., 51(7), 536581. https://doi.org/10.1016/10.1002/2015WR016921
 He and Sykes, [1996] He, Y., Sykes, J.F., 1996. On the spatialtemporal averaging method for modeling transport in porous media. Transport Porous Media, 22(1), 1–51. https://doi.org/10.1007/BF00974310
 Heße et al., [2009] Heße, F., Radu, F.A., Thullner, M. and Attinger, S., 2009. Upscaling of the advectiondiffusionreaction equation with Monod reaction. Adv. Water Resour., 32(8), 1336–1351. https://doi.org/10.1016/j.advwatres.2009.05.009
 Hornung, [1996] Hornung, U., 1996. Homogenization and porous media (Vol. 6). Springer Science & Business Media. https://link.springer.com/book/10.1007/9781461219200
 Illiano et al., [2020] Illiano, D., Pop, I.S., Radu, F.A., 2020. Iterative schemes for surfactant transport in porous media. Comput. Geosci. https://doi.org/10.1007/s10596020099492
 Irving and Kirkwood, [1950] Irving, J.H., Kirkwood, J.G., 1950. The statistical mechanical theory of transport processes. IV. The equations of hydrodynamics. J. Chem. Phys., 18(6), 817–829. https://doi.org/10.1063/1.1747782
 Kumar et al., [2011] Kumar, K., Van Noorden, T.L., Pop, I.S., 2011. Effective dispersion equations for reactive flows involving free boundaries at the microscale. Multiscale Model. Simul., 9(1), 29–58. https://doi.org/10.1137/100804553
 Kumar et al., [2016] Kumar, K., NeussRadu, M., Pop, I.S., 2016. Homogenization of a pore scale model for precipitation and dissolution in porous media. IMA J. Appl. Math., 81(5), 877–897. https://doi.org/10.1093/imamat/hxw039
 List and Radu, [2016] List, F. and Radu, F.A., 2016. A study on iterative methods for solving Richards’ equation. Comput. Geosci. 20 (2), 341–353. https://doi.org/10.1007/s1059601695663
 Mailloux et al., [2003] Mailloux, B.J., Fuller, M.E., Rose, G.F., Onstott, T.C., DeFlaun, M.F., Alvarez, E., Hemingway, C., Hallet, R.B., Phelps, T.J., Griffin, T., 2003. A modular injection system, multilevel sampler, and manifold for tracer tests. Ground Water, 41(6), 816–27. http://dx.doi.org/10.1111/j.17456584.2003.tb0242
 McClure et al., [2017] McClure, J.E., Dye, A.L., Miller, C.T., Gray, W.G., 2017. On the consistency of scale among experiments, theory, and simulation. Hydrol. Earth Syst. Sci., 21(2), 1063. 10.5194/hess2110632017
 Nolen et al., [2008] Nolen, J., Papanicolaou, G., Pironneau, O., 2008. A framework for adaptive multiscale methods for elliptic problems. Multiscale Model. Simul., 7, 171–196. http://dx.doi.org/10.1137/070693230
 Nordbotten et al., [2007] Nordbotten, J.M., Celia, M.A., Dahle, H.K., Hassanizadeh, S.M., 2007. Interpretation of macroscale variables in Darcy’s law. Water Resour. Res., 43(8), W08430. http://dx.doi.org/10.1029/2006WR005018
 Pop et al., [2004] Pop, I.S., Radu, F.A., Knabner, P., 2004. Mixed finite elements for the Richards’ equation: linearization procedure. J. Comput. Appl. Math., 168 (12), 365–373. https://doi.org/10.1016/j.cam.2003.04.008
 Porta et al., [2012] Porta, G.M., Riva, M., Guadagnini, A., 2012. Upscaling solute transport in porous media in the presence of an irreversible bimolecular reaction. Adv. Water Resour., 35, 151–162. https://doi.org/10.1016/j.advwatres.2011.09.004
 Schüler et al., [2016] Schüler, L., Suciu, N., Knabner, P., Attinger, S. 2016. A time dependent mixing model to close PDF equations for transport in heterogeneous aquifers. Adv. Water Resour., 96, 55–67. http://dx.doi.org/10.1016/j.advwatres.2016.06.012
 Schwarze et al., [2001] Schwarze, H., Jaekel, U., Vereecken, H., 2001. Estimation of macrodispersion by different approximation methods for flow and transport in randomly heterogeneous media. Transp. Porous Media, 43(2), 265287. https://doi.org/10.1023/A:1010771123844
 Schwede et al., [2008] Schwede, R.L., Cirpka, O.A., Nowak, W., Neuweiler, I., 2008. Impact of sampling volume on the probability density function of steady state concentration. Water Resour. Res., 44(12), W12433. http://dx.doi.org/10.1029/2007WR006668
 Skøien et al., [2003] Skøien, J.O., Blöschl, G., Western, A.W., 2003. Characteristic space scales and timescales in hydrology. Water Resour. Res., 39(10). http://dx.doi.org/10.1029/2002WR001736
 SoleMari and FernàndezGarcia, [2018] SoleMari, G., FernàndezGarcia, D., 2018, Lagrangian modeling of reactive transport in heterogeneous porous media with an automatic locally adaptive particle support volume, Water Resour. Res. 54(10), 8309–8331, http://dx.doi.org/10.1029/2018WR023033
 SoleMari et al., [2019] SoleMari, G., Bolster, D., FernàndezGarcia, D., SanchezVila, X., 2019. Particle density estimation with gridprojected and boundarycorrected adaptive kernels. Adv. Water Resour. 131, 103382. https://doi.org/10.1016/j.advwatres.2019.103382
 Srzic et al., [2013] Srzic, V., Cvetkovic, V., Andricevic, R., Gotovac, H., 2013. Impact of aquifer heterogeneity structure and localscale dispersion on solute concentration uncertainty. Water Resour. Res., 49(6), 3712–3728. http://dx.doi.org/10.1002/wrcr.20314
 Suciu, [2001] Suciu, N., 2001. On the connection between the microscopic and macoscopic modeling of the thermodynamic processes. Ed. Univ. Piteşti (in Romanian). English abstract: http://ictp.acad.ro/suciu/abstractPhDThesisSuciu.pdf.
 Suciu et al., [2015] Suciu, N., Radu, F.A., Attinger, S., Schüler, L., Knabner, P., 2015. A FokkerPlanck approach for probability distributions of species concentrations transported in heterogeneous media. J. Comput. Appl. Math., 289, 241–252. http://dx.doi.org/10.1016/j.cam.2015.01.030
 Suciu et al., [2016] Suciu, N., Schüler, L., Attinger, S., Knabner, P., 2016. Towards a filtered density function approach for reactive transport in groundwater. Adv. Water Resour., 90, 83–98. http://dx.doi.org/10.1016/j.advwatres.2016.02.016
 Suciu, [2019] Suciu, N., 2019. Diffusion in Random Fields. Applications to Transport in Groundwater. Birkhäuser, Cham. https://doi.org/10.1007/9783030150815
 Suciu and Radu, [2022] Suciu, N., Radu, F.A., 2022. Global random walk solvers for reactive transport and biodegradation processes in heterogeneous porous media. Adv. Water Resour. 166, 104268. http://dx.doi.org/10.1016/j.advwatres.2022.104268
 Suciu et al., [2021] Suciu, N., Illiano, D., Prechtel, A., Radu, F.A., 2021. Global random walk solvers for fully coupled flow and transport in saturated/unsaturated porous media. Adv. Water Resour., 152, 103935. https://doi.org/10.1016/j.advwatres.2021.103935
 [46] Suciu, N., Radu, F.A., Pop, I.S., 2023a. https://github.com/PMFlow/SpaceTimeUpscaling. Git repository.
 [47] Suciu, N., Radu, F.A., Pop, I.S., 2023b, Spacetime upscaling of reactive transport in porous media (extended version). Preprint arXiv.2112.10692. https://doi.org/10.48550/arXiv.2112.10692
 Teutsch et al., [2000] Teutsch, G., Ptak, T., Schwarz, R., Holder, T., 2000. Ein neues Verfahren zur Quantifizierung der Grundwasserimmission: I. Theoretische Grundlagen. Grundwasser, 4(5), 170–175. http://dx.doi.org/10.1007/s76700083687
 Tonina and Belin., [2008] Tonina, D., Bellin, A., 2008. Effects of porescale dispersion, degree of heterogeneity, sampling size, and source volume on the concentration moments of conservative solutes in heterogeneous formations. Adv. Water Resour., 31(2), 339–354. http://dx.doi.org/doi:10.1016/j.advwatres.2007.08.009
 Vamoş et al., [1996] Vamoş C., Georgescu A., Suciu N., Turcu I., 1996. Balance equations for physical systems with corpuscular structure. Phys. Stat. Mech. Appl., 227(1), 81–92. http://dx.doi.org/10.1016/03784371(95)003738.
 Vamoş et al., [1997] Vamoş C., Suciu N,, Georgescu, A., 1997. Hydrodynamic equations for onedimensional systems of inelastic particles. Phys. Rev. E, 55(5), 6277–6280. http://doi.org/10.1103/PhysRevE.55.6277
 Vamoş et al., [2000] Vamoş, C., Suciu, N., Blaj, W., 2000. Derivation of onedimensional hydrodynamic model for stock price evolution. Phys. Stat. Mech. Appl., 287(3), 461–467. http://dx.doi.org/10.1016/S03784371(00)00385X
 Vamoş et al., [2003] Vamoş, C., Suciu, N., Vereecken, H., 2003. Generalized random walk algorithm for the numerical modeling of complex diffusion processes. J. Comput. Phys., 186, 527–544. https://doi.org/10.1016/S00219991(03)000731
 Vamoş, [2007] Vamoş C., 2007. Continuous fields associated to corpuscular systems. Risoprint, ClujNapoca (in Romanian). English version: http://ictp.acad.ro/vamos/cvthesis/
 Van Duijn and van der Zee, [2018] Van Duijn, C.J., van der Zee, S.E.A.T.M., 2018. Large time behaviour of oscillatory nonlinear solute transport in porous media. Chemical Engineering Science, 183, 8694. https://doi.org/10.1016/j.ces.2018.02.045
 Wiedemeier et al., [1999] Wiedemeier, T.H., Rifai, H.S., Newell, C.J., Wilson, J.T., 1999. Natural attenuation of fuels and chlorinated solvents in the subsurface. John Wiley & Sons.
 Wright et al., [2021] Wright, E. E., Sund, N. L., Richter, D. H., Porta, G. M., Bolster, D., 2021. Upscaling bimolecular reactive transport in highly heterogeneous porous media with the LAgrangian Transport Eulerian Reaction Spatial (LATERS) Markov model. Stoch. Environ. Res. Risk. Assess., 35, 1529–1547. https://doi.org/10.1007/s0047702102006z
[1] Alecsa, C.D., Boros, I., Frank, F., Knabner, P., Nechita, M., Prechtel, A., Rupp, A., Suciu, N., 2019. Numerical benchmark study for fow in heterogeneous aquifers. Adv. Water Resour., 138, 103558. https://doi.org/10.1016/j.advwatres.2020.103558
[2] Andricevic R., 1998. Effects of local dispersion and sampling volume on the evolution of concentration fluctuations in aquifers. Water Resour. Res., 34(5):1115–29. http://dx.doi.org/10.1029/98WR00260
[3] de Barros, F.P.J., Fiori, A., 2014. Firstorder based cumulative distribution function for solute concentration in heterogeneous aquifers: Theoretical analysis and implications for human health risk assessment. Water Resour. Res., 50(5), 4018–4037. http://dx.doi.org/10.1002/2013WR015024
[4] Battiato, I. Multiscale dynamics of reactive fronts in heterogeneous media with fluctuating forcings. SIAM Conference on Mathematical and Computational Issues in the Geosciences September 11–14, 2017, Erlangen, Germany. https://archive.siam.org/meetings/gs17/program/invitedtalksprizepresentations/.
[5] BayerRaich, M., Jarsjo, J., Liedl, R., Ptak, T., Teutsch, G.2004. Average contaminant concentration and mass flow in aquifers from timedependent pumping well data: Analytical framework. Water Resour. Res., 40(8), W08303. http://dx.doi.org/10.1029/2004WR003095
[6] Bause, M., Knabner, P., 2004. Numerical simulation of contaminant biodegradation by higher order methods and adaptive time stepping. Comput. Visual. Sci., 7(2), 61–78. https://doi.org/10.1007/s007910040139y
[7] Berkowitz, B., 2021. HESS Opinions: Chemical transport modeling in subsurface hydrological systems – Space, time, and the holy grail of “upscaling”. Hydrol. Earth. Syst. Sci. Discuss., [preprint]. https://doi.org/10.5194/hess2021537
[8] Brunner, F., Radu, F.A., Bause, M., Knabner, P., 2012. Optimal order convergence of a modified BDM1 mixed finite element scheme for reactive transport in porous media. Adv. Water Resour., 35, 163–171. http://dx.doi.org/10.1016/j.advwatres.2011.10.001
[9] Bringedal, C., Berre, I., Pop, I.S., Radu, F.A., 2015. A model for nonisothermal flow and mineral precipitation and dissolution in a thin strip. J. Comput. Appl. Math., 289, 346–355. https://doi.org/10.1016/j.cam.2014.12.009
[10] Bringedal, C., Berre, I., Pop, I. S., Radu, F. A., 2016. Upscaling of nonisothermal reactive porous media flow under dominant Peclet number: the effect of changing porosity. Multiscale Model. Simul. 14(1), 502–533. https://doi.org/10.1137/15M1022781
[11] Caroni, E., Fiorotto, V., 2005. Analysis of concentration as sampled in natural aquifers. Transport Porous Media, 59(1),19–45. http://dx.doi.org/10.1007/s112420041119x
[12] Cirpka, O.A., Frind, E.O., Helmig, R., 1999. Numerical simulation of biodegradation controlled by transverse mixing. J. Contam. Hydrol., 40(2), 159182. https://doi.org/10.1016/S01697722(99)000443
[13] Cushman, J.H., 1983. Multiphase transport equations: Igeneral equation for macroscopic statistical, local spacetime homogeneity. Transport Theor. Stat. Phys., 12(1), 35–71. http://dx.doi.org/10.1080/00411458308212731
[14] Destouni, G., Graham, W., 1997. The influence of observation method on local concentration statistics in the subsurface. Water Resour. Res., 33(4), 663–676. http://dx.doi.org/10.1029/96WR03955
[15] Eberhard, J.P., Suciu, N., Vamos, C., 2007. On the selfaveraging of dispersion for transport in quasiperiodic random media. J. Phys. Math. Theor., 40(4), 597–610. https://doi.org/10.1088/17518113/40/4/002
[16] Gray, W.G., Miller, C.T., 2013. A generalization of averaging theorems for porous medium analysis. Adv. Water Resour. 62, 227–37. https://doi.org/10.1016/j.advwatres.2013.06.006
[17] Gray, W.G., Dye, A.L., McClure, J.E., PyrakNolte, L.J., Miller, C.T., 2015. On the dynamics and kinematics of twofluidphase flow in porous media. Water Resour. Res., 51(7), 536581. https://doi.org/10.1016/10.1002/2015WR016921
[18] He, Y., Sykes, J.F., 1996. On the spatialtemporal averaging method for modeling transport in porous media. Transport Porous Media, 22(1), 1–51. https://doi.org/10.1007/BF00974310
[19] Hornung, U., 1996. Homogenization and porous media (Vol. 6). Springer Science & Business Media. https://link.springer.com/book/10.1007/9781461219200
[20] Illiano, D., Pop, I.S., Radu, F.A., 2020. Iterative schemes for surfactant transport in porous media. Comput. Geosci. https://doi.org/10.1007/s10596020099492
[21] Klofkorn, R., Kr¨oner, D., Ohlberger, M., 2002. Local adaptive methods for convection dominated problems. Int. J. Numer. Meth. Fluid., 40(12), 79–91. https://doi.org/10.1002/fld.268
[22] Kumar, K., Van Noorden, T.L., Pop, I.S., 2011. Effective dispersion equations for reactive flows involving free boundaries at the microscale. Multiscale Model. Simul., 9(1), 29–58. https://doi.org/10.1137/100804553
[23] Kumar, K., NeussRadu, M., Pop, I.S., 2016. Homogenization of a pore scale model for precipitation and dissolution in porous media. IMA J. Appl. Math., 81(5), 877–897. https://doi.org/10.1093/imamat/hxw039
[24] List, F. and Radu, F.A., 2016. A study on iterative methods for solving Richards’ equation. Comput. Geosci. 20 (2), 341–353. https://doi.org/10.1007/s1059601695663
[25] Mailloux B.J., Fuller, M.E., Rose, G.F., Onstott, T.C., DeFlaun, M.F., Alvarez, E., Hemingway, C., Hallet, R.B., Phelps, T.J., Griffin, T., 2003. A modular injection system, multilevel sampler, and manifold for tracer tests. Ground Water, 41(6), 816–27. http://dx.doi.org/10.1111/j.17456584.2003.tb0242
[26] McClure, J.E., Dye, A.L., Miller, C.T., Gray, W.G., 2017. On the consistency of scale among experiments, theory, and simulation. Hydrol. Earth Syst. Sci., 21(2), https://doi.org/10.5194/hess2110632017
[27] Nolen, J., Papanicolaou, G., Pironneau, O., 2008. A framework for adaptive multiscale methods for elliptic problems. Multiscale Model. Simul., 7, 171–196. http://dx.doi.org/10.1137/070693230
[28] Nordbotten, J.M., Celia, M.A., Dahle, H.K., Hassanizadeh, S.M., 2007. Interpretation of macroscale variables in Darcy’s law. Water Resour. Res., 43(8), W08430. http://dx.doi.org/10.1029/2006WR005018
[29] Pop, I.S., Radu, F.A., Knabner, P., 2004. Mixed finite elements for the Richards’ equation: linearization procedure. Journal of computational and applied mathematics 168 (12), 365–373. https://doi.org/10.1016/j.cam.2003.04.008
[30] Porta, G.M., Riva, M., Guadagnini, A., 2012. Upscaling solute transport in porous media in the presence of an irreversible bimolecular reaction. Adv. Water Resour., 35, 151–162. https://doi.org/10.1016/j.advwatres.2011.09.004
[31] Schuler, L., Suciu, N., Knabner, P., Attinger, S. 2016. A time dependent mixing model to close PDF equations for transport in heterogeneous aquifers. Adv. Water Resour., 96:55–67. http://dx.doi.org/10.1016/j.advwatres.2016.06.012
[32] Schwarze, H., Jaekel, U., Vereecken, H., 2001. Estimation of macrodispersion by different approximation methods for flow and transport in randomly heterogeneous media. Transport Porous Media, 43(2), 265287. https://doi.org/10.1023/A:1010771123844
[33] Schwede, R.L., Cirpka, O.A., Nowak, W., Neuweiler, I., 2008. Impact of sampling volume on the probability density function of steady state concentration. Water Resour. Res., 44(12), W12433. http://dx.doi.org/10.1029/2007WR006668
[34] Skøien, J.O., Bloschl, G., Western, A.W., 2003. Characteristic space scales and timescales in hydrology. Water Resour. Res., 39(10). http://dx.doi.org/10.1029/2002WR001736
[35] Srzic, V., Cvetkovic, V., Andricevic, R., Gotovac, H., 2013. Impact of aquifer heterogeneity structure and localscale dispersion on solute concentration uncertainty. Water Resour. Res., 49(6), 3712–3728. http://dx.doi.org/10.1002/wrcr.20314
[36] Suciu N., 2001. On the connection between the microscopic and macoscopic modeling of the thermodynamic processes. Ed. Univ. Pitesti (in Romanian). English abstract: http://www1.am.unierlangen.de/suciu/Rabsth.pdf.
[37] Suciu, .N, Radu, F.A., Attinger, S., Sch¨uler, L., Knabner, P., 2015. A FokkerPlanck approach for probability distributions of species concentrations transported in heterogeneous media. J. Comput. Appl. Math., 289, 241–252. http://dx.doi.org/10.1016/j.cam.2015.01.030
[38] Suciu, N., Schuler, L., Attinger, S., Knabner, P., 2016. Towards a filtered density function approach for reactive transport in groundwater. Adv. Water Resour., 90, 83–98. http://dx.doi.org/10.1016/j.advwatres.2016.02.016
[39] Suciu, N., 2019. Diffusion in Random Fields. Applications to Transport in Groundwater. Birkhauser, Cham. https://doi.org/10.1007/9783030150815
[40] Suciu, N., Illiano, D., Prechtel, A., Radu, F.A., 2021. Global random walk solvers for fully coupled flow and transport in saturated/unsaturated porous media. Adv. Water Resour., 152, 103935. https://doi.org/10.1016/j.advwatres.2021.103935
[41] Suciu, N., Radu, F.A., 2021. Global random walk solvers for reactive transport and biodegradation processes in heterogeneous porous media. arXiv: 2107.08745
[42] Suciu, N., Radu, F.A., 2021. https://github.com/PMFlow/MonodReactions. Git repository. https://doi.org/10.5281/zenodo.5104076
[43] Suciu, N., Radu, F.A., 2021. https://github.com/PMFlow/SpaceTimeUpscaling. Git repository.
[44] Teutsch, G., Ptak, T., Schwarz, R., Holder, T., 2000. Ein neues Verfahren zur Quantifizierung der Grundwasserimmission: I. Theoretische Grundlagen. Grundwasser, 4(5):170–175. http://dx.doi.org/10.1007/s76700083687
[45] Tonina, D., Bellin, A., 2008. Effects of porescale dispersion, degree of heterogeneity, sampling size, and source volume on the concentration moments of conservative solutes in heterogeneous formations. Adv. Water Resour., 31(2), 339–354. http://dx.doi.org/doi:10.1016/j.advwatres.2007.08.009
[46] Vamos C., Georgescu A., Suciu N., Turcu I., 1996. Balance equations for physical systems with corpuscular structure. Phys. Stat. Mech. Appl., 227(1), 81–92. http://dx.doi.org/10.1016/03784371(95)003738
[47] Vamos C., Suciu N,, Georgescu, A., 1997. Hydrodynamic equations for onedimensional systems of inelastic particles. Phys. Rev. E, 55(5), 6277–6280. http://doi.org/10.1103/PhysRevE.55.6277
[48] Vamos, C., Suciu, N., Peculea, M., 1997. Numerical modelling of the onedimensional diffusion by random walkers. Rev. Anal. Numer. Theor. Approx., 26(1), 237–247. https://ictp.acad.ro/jnaat/journal/article/view/1997vol26nos12art32.
[49] Vamos, C., Suciu, N., Blaj, W., 2000. Derivation of onedimensional hydrodynamic model for stock price evolution. Phys. Stat. Mech. Appl., 287(3), 461–467. http://dx.doi.org/10.1016/S03784371(00)00385X
[50] Vamos, C., Suciu, N., Vereecken, H., 2003. Generalized random walk algorithm for the numerical modeling of complex diffusion processes. J. Comput. Phys., 186, 527–544. https://doi.org/10.1016/S00219991(03)000731
[51] Vamos C., 2007. Continuous fields associated to corpuscular systems. Risoprint, ClujNapoca (in Romanian). English version: https://ictp.acad.ro/vamos/cvthesis/
[52] Wiedemeier, T.H., Rifai, H.S., Newell, C.J., Wilson, J.T., 1999. Natural attenuation of fuels and chlorinated solvents in the subsurface. John Wiley & Sons.
[53] Wright, E. E., Sund, N. L., Richter, D. H., Porta, G. M., Bolster, D., 2021. Upscaling bimolecular reactive transport in highly heterogeneous porous media with the LAgrangian Transport Eulerian Reaction Spatial (LATERS) Markov model. Stoch. Environ. Res. Risk. Assess., 35, 1529–1547. https://doi.org/10.1007/s0047702102006z