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, Bayer-Raich 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 flux-weighted 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 non-smooth fields. This issue is overcome in kernel density estimations [Fernàndez-Garcia and Sanchez-Vila,, 2011, Sole-Mari and Fernàndez-Garcia,, 2018, Sole-Mari 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àndez-Garcia and Sanchez-Vila,, 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 space-time 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 [Bayer-Raich 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 [Bayer-Raich 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 space-time 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 well-posed 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 two-fluid 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 Navier-Stokes 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 spatio-temporal averaging, lead to closed upscaled equations with macroscopic coefficients completely determined by solving unit-cell 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 non-periodic and random parameters of the subsurface system [Nolen et al.,, 2008]. The concept of spatio-temporal upscaling in the context of homogenization by multiple-scale 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 spatio-temporal 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, Bayer-Raich et al.,, 2004] or to the spatial and temporal scales of the hydrological processes [Skøien et al.,, 2003].

In both direct spatio-temporal 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 spatio-temporal scales that can be directly related to those of the measurements and hydrological observations is the coarse-grained (CG) description of corpuscular systems through space-time (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 continuous-time random walk approach, with a power-law 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 spatio-temporal scales of interest. The fine-grained concentrations are represented microscopically by the actual number of molecules involved in reactions, their advective-diffusive movement being described by simple discrete-time random walk processes modeled with GRW algorithms [Suciu and Radu,, 2022]. This molecular description implicitly accounts for spatio-temporal 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 two-dimensional 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 space-time upscaling. Further, in Section 4 the CGST approach is applied to various one-dimensional problems of reactive transport solved by GRW algorithms, with emphasis on biodegradation processes in saturated and unsaturated natural porous media. Two-dimensional 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 space-time 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, ].

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 one-dimensional 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 $𝒩$ particles is described by piecewise analytic time functions ${\phi }_i(t):[0,T]⟼ℝ$, $i=1,\mathrm{\dots },𝒩$, then

1) there exists a macroscopic description of the same system given by almost everywhere (a.e.) continuous fields through coarse-grained space-time averages

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) $⟨\phi ⟩(𝐱,t;a,\tau )$ has continuous partial derivatives a.e. in ${ℝ}^d\times (\tau ,T-\tau )$ and satisfies the identity

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 }𝒩$, the CGST average (1) counts the number of particles inside the cube $C(𝐱,a)$ during a time interval of length $2\tau$, that is, $⟨1⟩$ is an a.e. continuous approximation of the macroscopic concentration field. The average over the number of particles defined by $\overline{𝝃}=⟨𝝃⟩/⟨1⟩$, if $⟨1⟩>0$, and by $\overline{𝝃}=0$, if $⟨1⟩=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 $𝐮$, according to (1), $\overline{𝝃}=𝐮$. In the absence of chemical reactions and jump discontinuities the term $\delta \phi$ in (2) vanishes and one obtains a continuity equation,

For $\tau ⟶0$ 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(𝐱,t)$ at the scale of the mathematical continuum, hereafter called “fine-grained concentration”, is defined as an average $M$ over the ensemble of realizations of the transport process,

where $\delta$ is the Dirac function and ${𝐗}_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

one obtains the following result.

The ensemble average of the coarse-grained average $⟨1⟩(𝐱,t;a,\tau )$ is an average of the fine-grained concentration $c(𝐱,t)$ over the cube $C(𝐱,a)$ and the time interval $(t-\tau ,t+\tau )$ given by

The ensemble average $M[⟨1⟩](𝐱,t;a,\tau )$ defines the continuous concentration field $c(𝐱,t;a,\tau )$ as observed at the spatio-temporal 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 ⟶0$, 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 ${𝐱}_i(t)$ and velocities ${𝝃}_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(𝐱,a)$, in the time interval $(t-\tau ,t+\tau )$.

The goal of the space-time 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 one-dimensional 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.

We consider the one-dimensional advection-diffusion transport in saturated porous media of two reactive chemical species governed at the fine-grained scale by the coupled system of equations

As shown in [Destouni and Graham,, 1997], the space-time 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 space-time upscaling by hypothetical experiments with independent spatial and temporal scales which ensure a good resolution of the averages and are bellow the space-time 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 $⟨1_{\nu }⟩$, at a given time point is quantified by relative errors computed with vector $l^2$-norms as $e_{c_{\nu }}=\Vert {\overline{c}}_{\nu }-⟨1_{\nu }⟩\Vert /\Vert ⟨1_{\nu }⟩\Vert$, as well as by maximum relative differences

The influence of $a$ and $\tau$ on discrepancy is analyzed for a specific problem in B.

The two-dimensional model for reactive transport with biodegradation reactions, similar to the one-dimensional model considered above, is governed by the equations

Equations (20-21) are parameterized with the van Genuchten-Mualem model (18-19) and the same Monod reaction system as in the previous sections. The two-dimensional CGST averaging procedure, first verified by computing the velocity components and the diffusion coefficients, will be used in the following to compute space-time upscaled concentrations.

The CGST-GRW approach proposed in this article performs the space-time upscaling adapted to the spatio-temporal 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 advective-diffusive 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 CGST-GRW 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 time-dependence 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 time-dependent flow velocity.

• 5.

  According to (4), the CGST average concentration is equivalent to a space-time average over the continuous field of the fine-grained 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.