Generalised Multi-Rate Models for conjugate transfer in heterogeneous materials

TL;DR

This paper introduces a generalized multi-rate model for conjugate heat and mass transfer in heterogeneous materials, capturing complex interactions without assuming local equilibrium or specific geometries, and enabling direct coefficient calculation.

Contribution

It develops a new GMRM framework that generalizes MRMT, allowing for arbitrary shapes and micro-scale computations, improving modeling flexibility and accuracy.

Findings

Derivation of a spectral decomposition-based GMRM model.

Demonstration that MRMT is a leading order approximation of GMRM.

Highlighting the limitations of simple re-scaling of transfer coefficients.

Abstract

We propose a novel macroscopic model for conjugate heat and mass transfer between a \emph{mobile region}, where advective transport is significant, and a set of \emph{immobile regions} where diffusive transport is dominant. Applying a spatial averaging operator to the microscopic equations, we obtain a \emph{multi-continuum} model, where an equation for the average concentration in the mobile region is coupled with a set of equations for the average concentrations in the immobile regions. Subsequently, by mean of a spectral decomposition, we derive a set of equations that can be viewed as a generalisation of the multi-rate mass transfer (MRMT) model, originally introduced by Haggerty & Gorelick. This new formulation does not require any assumption on local equilibrium or geometry. We then show that the MRMT can be obtained as the leading order approximation, when the mobile concentration is in local equilibrium. The new Generalised Multi-Rate Mode (GMRM) has the advantage of providing a direct method for calculating the model coefficients for immobile regions of arbitrary shapes, through the solution of appropriate micro-scale cell problems. An important finding is that a simple re-scaling or re-parametrisation of the transfer rate coefficient (and thus, the memory function) is not sufficient to account for the flow field in the mobile region and the resulting non-uniformity of the concentration at the interfaces between mobile and immobile regions.

Tables

Table 1. Table 1: Evaluation of λ i ​ k subscript 𝜆 𝑖 𝑘 \lambda_{ik} and β i ​ k subscript 𝛽 𝑖 𝑘 \beta_{ik} for simple geometries for which there is an analytical solution of the unsteady diffusion equation (see [ 47 ] for details). Here L 𝐿 L represents half the length of the layer (the domain in the layer goes from − L 𝐿 -L to L 𝐿 L ), R 𝑅 R is the radius of the sphere or cylinder and ζ k subscript 𝜁 𝑘 \zeta_{k} is the k 𝑘 k -th zero of the zero-order Bessel function of the first kind.

Geometry	${\lambda }_{ik}$	${\beta }_{ik}$
1d Layer	${\left(2k-1\right)}^2{\pi }^2\frac{{𝒟}_i}{4L_i^2}$	$\frac{8}{{\left(2k-1\right)}^2{\pi }^2}$
Cylinder	${\zeta }_k^2\frac{{𝒟}_i}{R^2}$	$\frac{4}{{\zeta }_k^2}$
Sphere	$k^2{\pi }^2\frac{{𝒟}_i}{R^2}$	$\frac{6}{k^2{\pi }^2}$

Equations

$$\text{Pe}=\frac{UL}{{\mathcal{D}}_m},$$

$${\text{Pe}}_i=\frac{U_i{R}_i}{{\mathcal{D}}_i}\ll 1,\forall i=1,\dots ,N_i,$$

$$\frac{\partial c_m}{\partial t}+\nabla \cdot \left(\mathbf{u}{c}_m-{\mathcal{D}}_m\nabla c_m\right)=0,\mathbf{x}\in {\Omega }_m.$$

$$\frac{\partial c_i}{\partial t}={\mathcal{D}}_i{\nabla }^2{c}_i,\mathbf{x}\in {\Omega }_i,i=1,\dots ,N_i.$$

$$c_i=c_m,{\mathcal{D}}_i\frac{\partial c_i}{\partial n}={\mathcal{D}}_m\frac{\partial c_m}{\partial n},\mathbf{x}\in \partial {\Omega }_i.$$

$$\overline{c_m}\left(\mathbf{x},t\right)=\underset{\Omega }{\int }{K}_{\Omega }(\mathbf{x}-{\mathbf{x}}^{\prime })c_m({\mathbf{x}}^{\prime },t)\text{d}V=K_{\Omega }\ast c_m$$

$$K_{\Omega }(\mathbf{x}-{\mathbf{x}}^{\prime })=\frac{1}{V}\{\begin{array}{ll}1& \forall {\mathbf{x}}^{\prime }\in \Omega \left(\mathbf{x}\right)\\ 0& \text{otherwise}\end{array}$$

$$\overline{c_m}\left(\mathbf{x},t\right)=K_{\Omega }\ast c_m=\frac{1}{V}\underset{\Omega }{\int }{c}_m({\mathbf{x}}^{\prime },t)\text{d}V$$

$$K_{{\Omega }_m}(\mathbf{x}-{\mathbf{x}}^{\prime })=\frac{1}{V_m}\{\begin{array}{ll}1& \forall {\mathbf{x}}^{\prime }\in {\Omega }_m\left(\mathbf{x}\right)\\ 0& \text{otherwise}\end{array}$$

$$\widetilde{c_m}\left(\mathbf{x},t\right)=K_{{\Omega }_m}\ast c_m=\frac{1}{V_m}\underset{{\Omega }_m}{\int }{c}_m({\mathbf{x}}^{\prime },t)\text{d}V$$

$$\overline{c_m}={\beta }_m\widetilde{c_m},$$

$${\beta }_m\frac{\partial \widetilde{c_m}}{\partial t}+\sum _{i=1}^{N_i}\frac{1}{V}\underset{\partial {\Omega }_i}{\int }{\mathcal{D}}_i\frac{\partial c_i}{\partial n}\text{d}S={\beta }_m\frac{\partial \widetilde{c_m}}{\partial t}+\sum _{i=1}^{N_i}{\dot{\mathcal{M}}}_i\left(t\right)=-\nabla \cdot {\mathbf{J}}_m,$$

$${\mathbf{J}}_m=\overline{\mathbf{u}{c}_m}-\overline{{\mathcal{D}}_m\nabla c_m},$$

$${\dot{\mathcal{M}}}_i\left(t\right)=\frac{1}{V}\underset{\partial {\Omega }_i}{\int }{\mathcal{D}}_i\frac{\partial c_i}{\partial n}\text{d}S.$$

$${\mathbf{J}}_{m,\text{eff}}={\mathbf{u}}_{\text{eff}}\widetilde{c_m}-{\mathcal{D}}_{m,\text{eff}}\nabla \widetilde{c_m},$$

$$\widetilde{c_i}\left(\mathbf{x},t\right)=\frac{1}{V_i}\underset{{\Omega }_i}{\int }{c}_i\left(\mathbf{x},t\right)\text{d}V,$$

$$\frac{\partial \widetilde{c_i}}{\partial t}=\frac{1}{V_i}\underset{\partial {\Omega }_i}{\int }{\mathcal{D}}_i\frac{\partial c_i}{\partial n}\text{d}S=\frac{\dot{\mathcal{M}}\left(t\right)}{{\beta }_i},$$

$${\beta }_m\frac{\partial \widetilde{c_m}}{\partial t}+\sum _{i=1}^{N_i}{\beta }_i\frac{\partial \widetilde{c_i}}{\partial t}=-\nabla \cdot {\mathbf{J}}_{m,\text{eff}}$$

$$c_i\left(\mathbf{x},t\right)={\psi }_i\left(\mathbf{x},t\right)+c_i^{\prime }\left(\mathbf{x},t\right)$$

$$\{\begin{array}{ll}{\nabla }^2{\psi }_i=0,& \mathbf{x}\in {\Omega }_i\\ {\psi }_i\left(\mathbf{x},t\right)=c_m\left(\mathbf{x},t\right),& \mathbf{x}\in \partial {\Omega }_i\end{array}$$

$$\{\begin{array}{ll}\frac{\partial c_i^{\prime }}{\partial t}-{\mathcal{D}}_i{\nabla }^2{c}_i^{\prime }=-\frac{\partial {\psi }_i}{\partial t},& \mathbf{x}\in {\Omega }_i\\ c_i^{\prime }=0,& \mathbf{x}\in \partial {\Omega }_i\end{array}$$

$$\underset{{\Omega }_i}{\int }{\nabla }^2{\psi }_{i}d{V}_i=\underset{{\Omega }_i}{\int }\nabla \cdot \nabla {\psi }_{i}d{V}_i=\underset{\partial {\Omega }_i}{\int }\frac{\partial {\psi }_i}{\partial n}\text{d}{S}_i=0,$$

$$c_i^{\prime }=\sum _{j=1}^{\infty }{c}_{ij}^{\prime }\left(\mathbf{x},t\right){\varphi }_{ij}\left(\mathbf{x}\right)$$

$$\{\begin{array}{ll}\mathcal{D}{\nabla }^2{\varphi }_{ij}={\lambda }_{ij}{\varphi }_{ij},& \mathbf{x}\in {\Omega }_i\\ {\varphi }_{ij}=0,& \mathbf{x}\in \partial {\Omega }_i\end{array}$$

$${\lambda }_{ij}^{\star }=\frac{{\lambda }_{ij}{L}_i^2}{{\mathcal{D}}_i}$$

$$\frac{\partial c_{ik}^{\prime }}{\partial t}={\lambda }_{ik}{c}_{ik}^{\prime }-\frac{1}{A_i}\frac{\partial }{\partial t}\underset{{\Omega }_i}{\int }{\psi }_i{\varphi }_{ik}\text{d}V$$

$$m_{ij}=\underset{{\Omega }_i}{\int }{\varphi }_{ij}\text{d}V,w_{ij}=\frac{m_{ij}}{A_i}$$

$$c_{ij}=\frac{1}{w_{ij}}{c}_{ij}^{\prime }+\frac{1}{m_{ij}}\underset{{\Omega }_i}{\int }{\psi }_i{\varphi }_{ij}\text{d}V$$

$$\frac{\partial c_{ik}}{\partial t}={\lambda }_{ik}\left(c_{ik}-{\psi }_{ik}\right)$$

$${\psi }_{ik}\left(t\right)=\frac{1}{m_{ik}}\underset{{\Omega }_i}{\int }{\psi }_i\left(\mathbf{x},t\right){\varphi }_{ik}\left(\mathbf{x}\right)\text{d}V$$

Full text

1 Introduction

Conjugate transfer in heterogeneous media is of pivotal importance for a wide range of applications ranging from dispersion of contaminants in aquifers [1, 2, 3, 4, 5] and stagnation/recirculation zones [6, 7, 8, 9] to heat transfer in granular media and suspension flows [10], or colloid interface reactions [11, 12]. In all these systems, we are faced with one (or more) flowing fluid exchanging mass or energy with a set of quiescent regions or impermeable inclusions, where diffusion can be assumed to be the dominant transport process. In this work we will refer to the first as mobile region and the latter as immobile regions. This terminology introduces a classification based on the mathematical modeling of regions rather than their physical meaning, and therefore allows to draw conclusions that are widely applicable to a class of problems. Similarly, we assume that heat and mass transfer processes obey the same governing equations (therefore we do not consider, for example, phase change or other critical phenomena).

While transport in weakly heterogeneous media can be accurately described using stochastic perturbative approaches [13] (see [14, 15, 16] for an extensive review), typical flow structures and exchange phenomena arising from strong heterogeneities (see for example [17, 18]) can not be captured by low order expansions. In fact, predictions from these methods show significant discrepancies when compared against observations from field experiments [19, 20], numerical simulations (for example [21]) and laboratory experiments [22].

To predict transport in strongly heterogeneous systems, a large number of methods have been developed, the most common of which are:

Integro-differential formulations [4] where the mass transfer to the immobile region is represented as the convolution of the concentration with an appropriate memory function over the past history of the system.

• 2.

The Multi-Rate mass transfer [1], which consist in modeling the transfer between mobile and immobile regions as a system of first order reactions.

• 3.

The continuous time random walk [13], where the movement of solute particles in the heterogeneous medium is represented as random walks in time and space.

Furthermore, it has been demonstrated that these methods are substantially equivalent [23, 13, 24] and a unified formulation based on the multi-rate mass transfer has been proposed [24]. This somewhat arbitrary choice was based on the sound basis that (i) the multi-rate mass transfer is generally more intuitive than the other methods and that (ii) it allows localisation. In the present work, we will add one further reason to motivate such choice: (iii) that the multi-rate mass transfer can be derived from the microscopic equations exactly, and intuitively interfaced with results from homogenisation (see [25, 26] for an extensive review of homogenisation theory).

However, accurate estimation for the closure parameters of the multi-rate mass transfer model is still a largely debated topic. Specifically, the multi-rate mass transfer model of Haggerty $\&$ Gorelick [1] requires a couple of parameters for each first order reaction:

$\alpha_{\text{HG}}$: the apparent exchange rate coefficient.

• 2.

$\beta_{\text{HG}}$: the capacity ratio.

It was suggested [1, 24] that while these parameters are indeed functions of other variables (like material and geometrical properties) at a fundamental level, they should really be considered as the fundamental coefficients for the model. A formal approach to obtain these coefficient consists in expressing the inter-region transfer as a memory term in the governing equations for the mobile region [4, 13]. Such term results from the convolution of the accumulation term with a memory function [11], which is then expanded in series of other functions (generally exponentials). The free parameters arising from this operation correspond to the parameters of Haggerty $\&$ Gorelick and they can be evaluated on the basis of analytical solutions for simple geometries [27]. However, one notorious limitation of such approach is the lack of theoretical basis to describe the dependence of the apparent exchange rate coefficients on the Reynolds number in the mobile region [28, 29]. In fact, several studies [30, 31, 32] showed that an exponential memory function is inadequate to describe the dependence on the flow rate. As a result, more complicated memory functions have been proposed as ad hoc solutions [33, 34, 35, 36], often based on the breakthrough curves and lacking any sort of physical connection with the underlying geometry or material properties. Therefore, calibration using laboratory experiments or numerical simulations [37, 6] and data fitting are often employed to obtain model parameters in practice. As a result, current mathematical formulations of multi-rate models still consider (at a macroscopic level) the concentration in the mobile region in equilibrium for what concern the inter-region exchange.

In this work we propose a novel general derivation of the multi-rate mass transfer model that address the following modeling issues:

• i

Providing a unique way of calculating the model parameters, like a set of equations that can be solved once for a whole class of problems.

• ii

Including the effect of advective transport on the conjugate transfer in a way that is mathematically formal and physically sound.

• iii

Derivation from first principles containing a limited and clear set of assumptions. This with the aim of facilitating any extension in future works.

This work is structured as follows. In section 2 we describe the microscopic equations and the approximations we employ. In section 3 we present the upscaling methodology in details and in Section section 4 we show how the model of Haggerty $\&$ Gorelick can be obtained as a zero-order approximation of our model. In section 4, we also present higher order models and we summarise the model parameters in section 6. We conclude in section 7 with an outlook to future extensions of the current model.

Additional details on the homogenisation procedures can be found in A.

2 Assumptions and microscopic equations

2.1 Heterogeneous domain

We consider the scenario presented in fig. 1. Let us consider a heterogeneous medium composed of a ”mobile” region and a number of ”immobile” zones. Therefore, $\Omega=\Omega_{m}\cup_{i=1}^{N_{i}}\Omega_{i}$, where $\Omega_{m}$ is the region occupied by the mobile region and $N_{i}$ is the number of inclusions. The mobile region is exchanging mass with the immobile regions through the inclusions’ boundaries $\partial\Omega_{i}$. We also assume that the regions $\Omega_{i}$ are completely included inside $\Omega$ and that they are disconnected (i.e., they only border on $\Omega_{m}$).

In the following we will assume that transport within inclusions $\Omega_{i}$ is dominated by diffusion, while on $\Omega_{m}$, advection might not be negligible. Thus, we can define a Peclet number:

[TABLE]

where $U$ is a characteristic velocity of the fluid, $L$ is a characteristic length and $\mathcal{D}$ is a diffusion coefficient. We will therefore assume that in the immobile regions:

[TABLE]

while no assumption is made on the Peclet number in the mobile region.

2.2 Microscopic governing equations

We assume that the concentration field $c_{m}\left(\bm{x},t\right)$ in the mobile region is obeying the advection-diffusion equation at the microscopic scale:

[TABLE]

Furthermore, we have $N_{i}$ diffusion equations for the concentrations in the immobile regions:

[TABLE]

We assume here the immobile diffusion coefficients $\mathcal{D}_{i}$ to be constant. This can be easily relaxed to smooth or piece-wise smooth coefficients and will be subject of future studies (by decomposing into coupled sub-regions).

At the interfaces $\partial\Omega_{i}$ we enforce continuity of fields and fluxes:

[TABLE]

This choice of boundary conditions implies that immobile regions do not exchange mass with each other, but they are only connected through the mobile region.

3 Upscaling methodology

3.1 Spatial filtering

Standard multi-continuum models [38] can be obtained from eq. 3 and 4, by applying a spatial filtering operator to the governing equation for the mobile region using a REV (Representative Elemetary Volume) $\Omega$ as support. We will assume that such REVs have a local periodic behaviour or, in other words, that their geometry changes very slowly with $\bm{x}$. Specifically, we assume that the number and geometrical configuration of the immobile regions included in a region $\Omega\left(\bm{x}\right)$ centred on $\bm{x}$ is essentially equivalent to that of a region $\Omega(\bm{x}+\delta\bm{x})$ for a sufficiently small $\delta\bm{x}$. This procedure produces fields that are much smoother than the original ones.

It is important to notice that in our model described by equations 3 and 4, the diffusive modes are mostly excited by the conjugate transfer with the immobile regions and not by source terms due, for example, to bulk reactions.

Furthermore, we can consider a macroscopic domain $\Omega_{\text{macro}}$ given by the union of a number of REVs $\Omega$. $\Omega_{\text{macro}}$ is taken sufficiently large to contain a large number of REVs, but sufficiently small to consider all those REVs as equivalent (i.e., disregarding any variation of the REVs geometry and material properties with $\bm{x}$). Therefore, it is possible to interchange between $\Omega$ and $\Omega_{\text{macro}}$ when computing averages without any loss of generality.

Thus, we define the volume average of $c_{m}$ in the region $\Omega\left(\bm{x}\right)$ (but it can be trivially extended to $\Omega_{\text{macro}}$) as the top-hat filter of volume $V=\int_{\Omega}\text{d}V$ centred on $\bm{x}$:

[TABLE]

Where $K_{\Omega}$ is a filtering Kernel and $*$ is the convolution operator. A typical Kernel that is widely used in fluid dynamics and multiphase flows is the top-hat [39, 40, 41]:

[TABLE]

Where $\Omega\left(\bm{x}\right)$ is a REV centred on $\bm{x}$. In this formulation both $\overline{c_{m}}$ and $c_{m}$ are both functions of the spatial coordinate $\bm{x}$. However, the integral operator results in $\overline{c_{m}}$ to be much smoother than $c_{m}$ and therefore we will consider as $\overline{c_{m}}$ does not depend on space at scales smaller than $\Omega$. Similarly, $V$ is also a slowly varying function of $\bm{x}$. We can therefore write an explicit expression for $\overline{c_{m}}$:

[TABLE]

As commonly done for multiphase and compressible flows, we also define the Favre top-hat Kernel [42] as the volume average over the mobile region $\Omega_{m}\left(\bm{x}\right)$ centred on $\bm{x}$ and of volume $V_{m}=\int_{\Omega_{m}}\text{d}V$:

[TABLE]

And therefore, the Favre averaged concentration $\widetilde{c_{m}}$ can be written as:

[TABLE]

The spatial filtering procedure is illustrated in fig. 2, where the resulting Favre averaged concentration is much smoother than the original one.

Generally, $\Omega$ is itself a function of the coordinate $\bm{x}$ but, since we assume symmetry of $\Omega$ under translation. the integral commutes with spatial derivatives and we will therefore omit its dependence of $\bm{x}$ for the sake of brevity. such that we obtain the relation:

[TABLE]

where we introduced the capacity of the mobile region $\beta_{m}$ defined as $V_{m}/V$. Some authors (for example [1]) define $\beta_{m}$ as being multiplied by a retardation factor obtained from a re-scaling of the time coordinate. Without loss of generality, we will not consider the retardation factor explicitly.

Thus, the main difference between volume averaging and Favre averaging is that the first is carried over all the space (mobile and immobile regions), while the second is restricted to a particular region. The main reason to introduce such difference, is that it formally leads to equation 11, and thus to the definition of capacity.

3.2 Multi-continuum formulation

Assuming that the the immobile regions $\Omega_{i}$ are fully included into $\Omega$ (i.e., $\partial\Omega\cap\partial\Omega_{i}=\emptyset\;i=1,\dots,N_{i}$), applying the integral operator (6) to (3) and making use of the Green’s theorem we obtain:

[TABLE]

where we have defined the total average flux in the mobile region:

[TABLE]

and the average inter-region mass exchange rate for region $i$:

[TABLE]

Since in this work our focus is on the interface exchange, we do not perform an accurate upscaling of $\mathbf{J}_{m}$, and we will use a simplified expression without any loss of generality:

[TABLE]

where $\bm{u}_{\text{eff}}$ and $\mathcal{D}_{m,\text{eff}}$ are the effective velocity and the effective diffusivity in the mobile region. The capacity $\beta_{m}$ does not appear explicitly into (15) since it is generally accounted for within the effective parameters.

We then define the Favre averaged concentration in the immobile regions as:

[TABLE]

where $V_{i}=\int_{\Omega_{i}}\text{d}V$ is the volume occupied by region $\Omega_{i}$. Thus, we integrate (4) to obtain:

[TABLE]

where $\beta_{i}=V_{i}/V$ is the capacity of immobile region $i$. The time derivative in eq. 17 is a partial derivative since $c_{i}$ depends on $\bm{x}$ at the macro scale (for example, due to the distribution of immobile regions at the macroscale). eq. 17, substituted into (12), leads to the multi-continuum equation for the concentration field in the mobile region:

[TABLE]

In (18), we transformed the boundary conditions of the microscopic equation into source terms, one for each immobile region. However, in this formulation $\widetilde{c_{i}}$ still needs to be found through an equation that is valid at the microscopic scale and, thus, requires the complete knowledge of the concentration in the immobile region.

3.3 Multi-rate mass transfer

In order to express $\widetilde{c_{i}}$ in a closed form that only depends on the geometrical and physical properties of the immobile region (as well as from the boundary value of $c_{m}$), we perform the following decomposition:

[TABLE]

Where the function $\psi_{i}$ satisfies the following equation and boundary conditions:

[TABLE]

While $c_{i}^{\prime}$ is given by:

[TABLE]

Summing eq. 20 and 21, and using decomposition 19 gives back eq. 4 with the correct boundary conditions. Due to eq. 20, $\psi_{i}$ satisfies the following Gauß-Green integral:

[TABLE]

being $n$ a field normal to $\partial\Omega_{i}$ and $S_{i}$ the surface of $\partial\Omega_{i}$.

Therefore, in our formulation the function $\psi_{i}$ is simply required to satisfy the non-homogeneous time dependent boundary conditions, while $c_{i}$ satisfies a non-homogeneous unsteady diffusion equation with homogeneous boundary conditions.

The homogeneous form of eq. 21 leads to an eigenvalue problem following a separation of variables, and can therefore be expressed in series of eigenfunctions without any loss of generality:

[TABLE]

Where $c_{ij}^{\prime}\left(t\right)$ are series coefficients that depend on time and on $\bm{x}$ at the macro scale only, while the eigenfunctions $\phi_{ij}\left(\bm{x}\right)$ carry the dependence on the spatial coordinate at the micro scale and satisfy:

[TABLE]

Where $\lambda_{ij}$ is the eigenvalue corresponding to eigenfunction $\phi_{ij}$.

While our decomposition of the spatial dependence in 23 may look arbitrary at first sight, in practice it simply mean that there co-exist two problems for the immmobile regions: (i) a local one and (ii) a global one. The local problem (i) refers to the solution within the single immobile regions and is described by eq. 4 within the REV $\Omega$. The global problem (ii) involves how the fields in the immobile regions vary at a macroscopic scale and how the immobile regions communicate. In our case, the immobile regions are disconnected and therefore, the spatial dependence of $c_{ij}^{\prime}$ is only keeping track of the different initial conditions at the macroscopic scale (since the boundary conditions are homogeneous). Therefore, we will take $c_{ij}^{\prime}$ out of any spatial derivative within the immobile regions, since they are assumed to be negligible.

Both eq. 20 and 24 can be made dimensionless by rescaling with respect to a characteristic length of the inclusion $L_{i}$ and the diffusion coefficient in the immobile region $\mathcal{D}_{i}$. As a consequence, we can relate the dimensional eigenvalue $\lambda_{ij}$ with a dimensionless eigenvalue:

[TABLE]

Following this rescaling, eq. 20 and 24 are not just valid for a particular geometry, but for class of similar geometries.

Substituting solution 23 back into eq. 21 and projecting into $\phi_{ik}$ we obtain:

[TABLE]

Where $A_{i}=\int_{\Omega_{i}}\phi_{ik}\phi_{ik}\text{d}V$ is the normalisation factor of the eigenproblem, which depends on the geometry only.

For reasons that will be clear in the next section, we introduce the following definitions:

[TABLE]

[TABLE]

Substituting into eq. 26 we obtain:

[TABLE]

Where we introduced the projection of $\psi$ into $\phi_{ik}$ scaled over the norm of $\phi_{ik}$:

[TABLE]

Therefore, $c_{i}$ is now given by:

[TABLE]

Where:

[TABLE]

is the correction function for the immobile region, which accounts for the non-homogeneity of $c_{m}\left(\bm{x},t\right)$ at the interface.

3.3.1 Computation of the exchange rate

We now compute the Favre averaged concentration in the $i$-immobile region:

[TABLE]

Where the favre averaged correction function is given by:

[TABLE]

The terms $\beta_{ik}$ play the role of capacities (or a normalised weighting function) since:

[TABLE]

Equation 35 is the so-called partition of unity (notice that $\beta_{ik}$ is generally still functions of the spatial coordinates at the macroscale) and it arises directly from the eigenproblem. Recalling eq. 17, we then obtain an expression for the mobile-immobile exchange rate:

[TABLE]

It is important to notice that all the terms involved in the multi-rate transfer can be computed a priori by solving a cell problem, which consists in solving eq. 20 for $\psi_{i}$ and the eigenvalue problem 24 for each immobile region $i$. However, $\psi_{i}\left(\bm{x},t\right)$ is a non trivial function of $c_{m}\left(\bm{x},t\right)$, and in the present formulation its computation requires the solution of equation 20 for each instant of time. This is clearly not desirable, since it would mean that a numerical algorithm would have to solve eq. 20 at each time step. Furthermore, no information regarding the functional dependence of $c_{m}$ on the flow rate is provided in the current formulation. Therefore, we need to introduce some information regarding the micro-structure of $c_{m}\left(\bm{x},t\right)$ in order to make any further progress.

3.4 Representation of $c_{m}\left(\bm{x},t\right)$ using homogenisation theory

So far, our formulation is exact, in the sense that we made no assumption regarding the regularity of the fields and we retained all the terms arising from the volume averaging. However, we still do not have an expression for the concentration field at the interface between mobile and immobile regions since that would require the complete knowledge of $c_{m}\left(\bm{x},t\right)$.

In order to give a representation of the spatial variability of $c_{m}\left(\bm{x},t\right)$ without having to solve the microscopic unsteady governing equations, one can employ the classical two-scale expansion of homogenisation theory [26, 43], and express $c_{m}$ as:

[TABLE]

where $\bm{\chi}_{n}$ is the corrector function corresponding to the $n$-order of the series and $\mathrel{\mathop{\mathchar 58\relax}}$ represents the contraction between the corrector and the $n$-order gradient $\bm{\nabla}^{n}$.

$\widetilde{c_{m}}$ and $\bm{\nabla}^{n}\widetilde{c_{m}}$ varies much slower than $\bm{\chi_{n}}$ in $\bm{x}$ and can be considered as constant111They are, in fact, constant at the microscale, when a two-scale hypothesis is introduced. when plugged into the microscopic equations. One important feature of corrector tensors is that they provide crucial information on the transport anisotropy. For example, if the flow field is unidirectional, the first order corrector $\bm{\chi}_{1}$ will be a vector field oriented towards the flow direction but (unlike the velocity field) it will not be zero at the interface between mobile and immobile regions. Therefore, employing eq. 37 allows to reconstruct the interface concentration from the gradients of $\widetilde{c}_{m}$ weighted with functions of the transport properties, provided that expansion 37 is shown to be convergent (which is beyond the scope of this work). We illustrate how the correctors and the above expansion can be obtained using homogenisation theory in A.

As a side note, we mention that homogenisation can be also employed to obtain an expression for $\mathbf{J}_{m,\text{eff}}$ [44, 45] and it is therefore synergic to the current problem. Furthermore, homogenisation theory can also be employed in place of volume averaging to derive dual porosity models [46].

To introduce the information provided by the corrector equation into our problem, we can expand $\psi_{i}$ in a similar fashion:

[TABLE]

where the functions $\bm{\Psi_{n}}$ are coupled with the correctors $\bm{\chi}_{n}$ at the interface and satisfy (owing the linearity of equation 20) :

[TABLE]

These are a set of partial differential equations for tensors of rank $n$. Notice that $\bm{\Psi}_{in}$ satisfies a boundary integral relation similar to eq. 22. We can now substitute expansion 38 into $\widetilde{\theta_{i}}$ to obtain:

[TABLE]

where we introduced $\bm{\Theta}^{\star}_{in}$ as the internal corrector tensor of rank $n$ for immobile region $i$, which accounts for the internal effects of the spatial variability of $c_{m}\left(\bm{x},t\right)$ at the interface. This formulation shows that, when we assume $c_{m}\left(\bm{x},t\right)=\widetilde{c_{m}}$ at the interface, then $\widetilde{\theta_{i}}=0$, which means that no correction is necessary.

Expansion 38 can now be substituted into the evolution equation for $c_{ik}$, leading to:

[TABLE]

4 Governing equations of the generalised multi-rate transfer model

We can now write down a set of equations for the Generalised Multi-Rate Transfer (GMRT) model which can be closed using a set of parameters corresponding to different geometries. When a specific geometry is selected, such parameters are constants or are simple function of geometrical and material properties through a rescaling (as for the eigenvalue $\lambda_{ik}$ =$\lambda_{ik}^{\star}\mathcal{D}_{i}/L_{i}^{2}$). The full system of equations is:

[TABLE]

In this system (GMRT-TS) mixed time-space derivatives are present. In the next section, these will be replaced to obtain a more convenient form. Nevertheless, GMRT-TS is exact as long as $c_{m}\left(\bm{x},t\right)$ can be expanded using the corrector eq. 37, and the series is convergent. In practical applications, one would also truncate both the series in $n$ and $k$ to achieve the desired accuracy or retain only a certain number of terms. In that case, some considerations on the approach of the series to convergence are required. However, if the macroscopic field $\widetilde{c}_{m}$ is sufficiently regular it is possible to obtain a good approximation just with the first order corrector $\bm{\chi}_{1}$ [26].

A key feature of the current formulation is that accounts for the non-uniform distribution of the concentration field at the interface from a microscopic perspective and shows how this can be upscaled to a macroscopic set of equations. Surprisingly, this does not lead to a new exchange rate (which is equivalent to the eigenvalue of the homogeneous problem $\lambda_{ik}$), but instead to an additional term in the equations for $c_{ik}$ and a new rate term. These terms lead to mixed and potentially high order derivatives in the governing equation for the mobile concentration. However, in practical applications one rarely goes beyond a second order corrector and therefore this does not alter the order of the differential equation.

4.1 A note on the truncation of the multi-rate series: equilibrium modes

Clearly, practical applications require the multi-rate series to be truncated at some value $k_{\text{max}}$ corresponding to $\beta_{ik_{\text{max}}}$, $\lambda_{ik_{\text{max}}}$, $c_{ik_{\text{max}}}$, and $\bm{\Psi}_{ik_{\text{max}}n}$. While previous works enforced $\sum_{k}\beta_{ik}=1$ by rescaling the capacities [24], here we propose a more accurate and rigorous approach based the above mathematical derivation of the GMRT.

It is easy to see that the truncated modes $k>k_{\text{max}}$ will approach equilibrium faster, since they correspond to small perturbations inside the immobile region and to larger eigenvalues. They can be therefore assumed to be in equilibrium if $k_{\text{max}}$ is chosen sufficiently large, according to the physics of the problem. Therefore, one can write:

[TABLE]

Furthermore, the corresponding eigenfunctions $\phi_{ik}$ will be highly oscillating for $k>k_{\text{max}}$ in contrast with the corrector $\bm{\Psi}_{in}$, which we choose to be a smooth function (recall equation 20). Therefore the projection of $\bm{\Psi}_{in}$ over $\phi_{ik}$ will be very small for $k>k_{\text{max}}$ since $\bm{\Psi}_{in}$ will not have significantly high modes. The regularity of $\bm{\Psi}_{in}$ is also connected to the change of $c_{m}\left(\bm{x},t\right)$ over $\partial\Omega_{i}$. In most of the applications (e.g., forced convection) $c_{m}\left(\bm{x},t\right)$ varies regularly over the interfaces, and thus the correctors $\bm{\chi}_{n}$ vary smoothly (and slowly) over $\partial\Omega_{i}$. Therefore, we can also assume that there exist a $k_{\text{max}}$ such that:

[TABLE]

and, consequently,

[TABLE]

which means that for sufficiently large $k$ the modes are in equilibrium with the average field. Therefore, a more sensible approximation of the truncated terms would be a scaling of $\beta_{i}$ (and thus $\beta_{m}$) to account for the removed modes compared to simply rescaling $\beta_{ik}$. In practice, one would then have new truncated immobile capacities $\beta_{i}^{\text{tc}}$ and mobile capacities $\beta_{m}^{\text{tc}}$ defined as:

[TABLE]

Thus, the system behaves as if the mobile region was larger and the immobile regions were smaller (in terms of volumetric occupation, not geometrical parameters). This is simply due to the fact that we assumed that the dynamics of modes $k>k_{\text{max}}$ in the immobile regions is completely determined by the mobile region.

4.2 The multi-rate model of Haggerty & Gorelick: the leading order approximation

The original Multi-Rate Mass Transfer (MRMT) model proposed by Haggerty & Gorelick [1] can be obtained as a special case of our general formulation. More specifically, their model can be considered as a leading order approximation for $\widetilde{\psi_{i}}$, which results in the system:

[TABLE]

Therefore, the model of Haggerty & Gorelick is obtained under the approximation that the concentration at the interface between each immobile region and the mobile region is uniform and equal to $\widetilde{c_{m}}$. This is acceptable for systems where the mobile region is approximatively in local equilibrium at the microscale. This can be the case of a well-mixed concentration in the mobile region.

4.3 Computation of $\beta_{ik}$ and $\lambda_{ik}$

Coefficients $\beta_{ik}$ and $\lambda_{ik}$ do not depend in any way on the interface concentration $c_{m}\left(\bm{x},t\right)$ and, following our approach, they bear no dependence on the transport processes happening in the mobile region. Therefore, they can be calculated exactly using only geometrical shape and material properties of the immobile regions as input.

Table 1 shows the expression of $\lambda_{ik}$ and $\beta_{ik}$ for a set of simple geometries. Clearly, our coefficients match those proposed by Haggerty $\&$ Gorelick [1], except for a factor $\beta_{m}$ in $\beta_{ik}$, which is consistent with our formulation since we do not divide the equation for $\widetilde{c}_{m}$ by $\beta_{m}$.

In this approximation. coefficients $\lambda_{ik}$ and $\beta_{ik}$ have the same meaning as in Haggerty $\&$ Gorelick, where $\lambda_{ik}$ plays the role of exchange rate between $c_{ik}$ and $\psi_{ik}$. As demonstrated in table 1, $\lambda_{ik}$ is a function of geometrical dimensions and material properties through $\lambda_{ik}$ =$\lambda_{ik}^{\star}\mathcal{D}_{i}/L_{i}^{2}$, where the dimensionless eigenvalue depends on the shape of the immobile region only. On the contrary, $\beta_{ik}$ is a dimensionless weight that depends only on the class of geometrical shapes.

5 Beyond classic MRMT

While eq. 42 allows to easily recover the standard MRMT model in the limit of equilibrium concentration in the mobile region, the presence of a mixed derivative makes its physical interpretation rather cumbersome. Furthermore, such term can introduce instabilities in numerical solution algorithms.

To this end, it is useful to rewrite eq. 36 using the integral form of the exchange rate:

[TABLE]

Integrating the eigenvalue eq. 24 over $\Omega_{i}$, we can obtain the following relation for the eigenvalues:

[TABLE]

Expanding the first term on the right hand side of eq. 48 leads to:

[TABLE]

where we employed eq. 22 on the right-hand-side. Then, substituting expansion eq. 38 results into:

[TABLE]

The additional terms are consistent with the evolution equation for $c_{ik}$, so that the mobile-to-immobile fluxes are identical to the corresponding immobile-to-mobile flux regardless the number of terms retained in the expansions.

5.1 Generalised Multi-Rate Transfer equations

The complete set of equations 42 can be therefore rewritten without mixed time-space derivatives as:

[TABLE]

where we defined the equilibrium concentration for term $k$ of region $i$ as:

[TABLE]

This system of equations does not pose any significant issue for corrections up to the second order, since the order of the differential operators remains unchanged and no mixed derivatives arise. Physically, these additional terms change the equilibrium concentration at which $\partial c_{ik}/\partial t=0$.

5.2 First order correction and drift flux approximation

Retaining first order corrections in eq. 52 is equivalent to adding a drift like term to the standard multi-rate equation for the mobile region. The governing equations are given by:

[TABLE]

For the special case in which the material microstructure does not vary in space and the flow field is macroscopically homogeneous (i.e., $\bm{u}_{\text{eff}}=\text{const}$), $\bm{\Psi}_{ik1}$ does not depends on the spatial coordinates and we can define a drift velocity:

[TABLE]

and thus a new effective velocity:

[TABLE]

Therefore, the equation for the mobile region simply reduces to a standard advection diffusion equation, with an additional multi-rate reactive term:

[TABLE]

5.3 Physical considerations on $\bm{\Psi}_{ik1}$

eq. 52 is describing a reactive system where the equilibrium concentrations of the immobile regions are not the same as the concentration in the mobile region. Thus, in our model, the equilibrium point is shifted by the correctors, based on the gradients of $c_{m}$. This is a direct consequence of non-equilibrium at the microscale (i.e., within $\Omega$) and can be attributed (at least asymptotically) to the the flow field and to the existence of boundary layers, which effectively results in a different equilibrium concentration for each $c_{ik}$. Our approach based on the synergy between homogenisation theory and spectral decomposition provides a formal way to account for such non-equilibrium.

In order to understand the meaning of these corrector terms, it is useful to consider the toy case depicted in fig. 3.

Such system is fundamentally monodimensional, and can be characterised by having:

[TABLE]

Now, we consider eq. 37 at the first order:

[TABLE]

Considering the local concentration in the mobile region, if all immobile regions have the same initial concentration, it follows that the local maxima will be locate at the interfaces as depicted in fig. 3. Thus, considering that inequality 58 holds, the same inequality holds for the $x$ component of the first order corrector $\chi_{1,x}$:

[TABLE]

Therefore, if all components of $\bm{\nabla}\widetilde{c}_{m}$ are positive, all the components of $\bm{\chi}_{1}$ are also positive. This positivity is transferred to $\bm{\Psi}_{i1}$ through eq. 39 due to the properties of elliptic operators. While $\bm{\Psi}_{ik1}$ is not necessarily positive, the projection on the first eigenfunction $\bm{\Psi}_{i11}$ is positive.

It is easy to demonstrate that this positivity property holds also when $\bm{\nabla}\widetilde{c}_{m}<0$.

Therefore, as illustrated in fig. 3, a in a system with $\bm{\nabla}\widetilde{c}_{m}>0$, the equilibrium concentration in the immobile regions will be larger than $\widetilde{c}_{m}$ due to the higher value of $c_{m}$ at the interface. On the contrary, in the case $\bm{\nabla}\widetilde{c}_{m}<0$, this will be lower.

While such argument was based on the analysis of a simple system, it is often valid for a large range of situations as, for example, in aquifer remediation and in many applications it is possible to guess the sign of the correctors by looking at the gradients.

However, when different immobile regions have different initial conditions or the transfer in the immobile regions is strongly asymmetric, this positivity condition may be violated.

5.4 Second order correction and diffusive flux approximation

We now consider correction terms up to second order. Such term brings a second order differential operator into eq. 52:

[TABLE]

Now, we can decompose tensor $\bm{\Psi}_{ik2}$ into hydrostatic and deviatoric components:

[TABLE]

Where $\text{tr}\left(\bm{\Psi}_{ik2}\right)$ is the trace of $\bm{\Psi}_{ik2}$ and $\mathbf{I}$ is the identity tensor. We now introduce the diffusion coefficient arising from the conjugate transfer $\mathcal{D}_{\text{ct}}$:

[TABLE]

which correspond to the second order correction arising in the case of macroscopically isotropic material with istropic immobile regions.

Again, we make the approximation of homogeneous, isotropic material with macroscopically homogeneous velocity field so that $\mathcal{D}_{\text{ct}}$ does not depend on the spatial coordinate. Under these approximations, the second order correction term becomes a purely diffusive contribution and we can thus define a new total diffusion coefficient:

[TABLE]

Therefore, the equation for the mobile concentration simplifies to:

[TABLE]

6 Summary of model parameters

Clearly, the the multi-rate series would be generally truncated at a desired accuracy. All the correction terms arising in the formulation of the present model can be evaluated based on analytical or numerical analysis of the immobile and mobile regions. All such parameters can be evaluate a priori and do not require additional computation when solving the macroscopic problem.

Specifically, two parameters are independent on the flow and geometry in the mobile region:

• $\lambda_{ik}$

: these are simply the eigenvalues corresponding to the homogeneous eigenproblem in the immobile region.

• $\beta_{ik}$

: these weights can be calculated similarly to $\lambda_{ik}$, from the solution of the eigenproblem. Once the eigenfunctions are known, $\beta_{ik}$ is given by: $\beta_{ik}=(\int_{\Omega_{i}}\phi_{ik}\text{d}V)^{2}/(V_{i}\int_{\Omega_{i}}\phi_{ik}^{2}\text{d}V)$

These are the same parameters of standard multi-rate models. Furthermore, there ore other parameters that require the solution of a cell problem in the mobile region and therefore, that bring information regarding the interplay of conjugate transfer and transport in the mobile region. Such terms make use of the correctors $\bm{\chi}_{in}$ obtained from homogenisation theory.

• $\bm{\Psi}_{ikn}$

: Projection of the function $\bm{\Psi}_{in}$ on the eigenfunction $\phi_{ik}$ scaled with the norm of $\phi_{ik}$. Clearly, the number of these parameters equals the number of terms in the multi-rate expansion but one can exploit some knowledge of the microstructure to simplify their expression.

Other quantities we introduced, like $\mathcal{D}_{\text{ct}}$, $c^{\text{eq}}_{ik}$ or $\bm{u}_{\text{drift}}$, can be obtained from the other parameters.

It is worth to notice that, as it is often suggested for the standard MRMT, it is possible to consider each of the parameter as unknown and obtainable (for example) trough inverse analysis or data fitting. In this case, while the details of the derivation of $\lambda_{ik}$, $\beta_{ik}$ and $\bm{\Psi}_{ikn}$ become irrelevant, it is still crucial to remember that all the physics of non equilibrium is contained in $\bm{\Psi}_{ikn}$.

7 Conclusions

In this paper we propose a novel approach to derive the multi-rate mass transfer model that is different from that of the memory function or that of Haggerty $\&$ Gorelick. Our model is derived starting from the microscopic equations and it is parameter free, i.e., it is possible to directly evaluate all the closure parameters in a unique manner. While our method agrees with previous results obtained by Haggerty $\&$ Gorelick, it also contains their multi-rate model as a special case and allows extension to non-equilibrium situations, where the concentration in the mobile region is not uniform. Especially, when homogenisation techniques are employed to evaluate the effective transport in the mobile region, our method provides an exact framework for the upscaling of the conjugate transfer problem, the accuracy of which is given by the terms retained from the infinite series.

Our model predicts that additional arise in the governing equations of the multi-rate mass transfer when accounting for the effect of transport processes in the mobile region on the inter-region exchange. These terms are brought into the framework by the corrector equation resulting from homogenisation, which at the second order have the form of a drift and a diffusive contribution.

Furthermore, under the assumptions of isotropy and homogeneity these terms can be absorbed into the effective diffusivity and effective velocity, thus leaving the form of the governing equations in the mobile region unchanged. However, the concentration in each immobile region will now depend on high order spatial derivatives of the concentration in the mobile region.

Despite the self-consistency of this model (all the parameters can be evaluated from first principles without calibration) and its completeness with respect to the initial hypothesis (we never introduced additional hypothesis or simplifications in the development of our formulation) there are still some significant phenomena that should be accounted for when modelling real systems. Some examples are:

Exchange between immobile regions.

• 2.

Multiple mobile regions with different mobility (e.g., fractures).

• 3.

Chemical reactions at interfaces.

• 4.

Multiphase flow, heat and mass transfer.

• 5.

Internal flow currents in the immobile regions.

Future works could focus on one or more of these topics to improve the range of applicability of this proposed model.

8 Acknowledgements

This work has been funded by the European Union’s Horizon 2020 research and innovation programme, grant agreement number 764531, ”SECURe – Subsurface Evaluation of Carbon capture and storage and Unconventional risks”.

Appendix A Homogenisation: evaluation of the first order corrector

Homogenisation is a perturbative method that allows to separate the original multiscale problem into a hierarchy of problems acting at different scales. In the following, we show how an equation for the immobile concentration similar to equation 12 can be obtained using homogenisation and how to calculate the first order corrector $\bm{\chi}_{1}$. The purpose of this Appendix is simply to illustrate the method applied to the current study. For a detailed and rigorous description of the homogenisation procedure see for example [45].

The process starts defining an expansion parameter:

[TABLE]

and a microscopic scale $\bm{y}$ such that:

[TABLE]

Where $R$ is a characteristic length of the immobile regions and $L$ is a characteristic length at the macroscale. The field $c_{m}\left(\bm{x},\bm{y},t\right)$ is then expanded in asymptotic series of $\varepsilon$:

[TABLE]

Spatial differential operators are expanded to account for the microscopic scale:

[TABLE]

Here $\bm{y}$ represents the variation across the REV $\Omega$, while $\bm{x}$ is a coordinate on the macroscopic volume $\Omega_{\text{macro}}$. This splitting is also known as two-scale asymptotics.

Then, 68 and 69 are substituted into equation 12 and, retaining terms up to $\mathcal{O}\left(\varepsilon^{-2}\right)$:

[TABLE]

The boundary conditions are of the second type, in agreement with equation 5:

[TABLE]

Where the order $\varepsilon$ is taken due to the scaling of $\bm{f}_{i}$ with the specific surface and the negative sign takes into account for the vectors normal to the surface. This can be considered as an approximation of ”slow flux”, we are fundamentally assuming that the diffusion is dominant at the macroscale with respect to the inter-region flux.

The boundary conditions are also expanded (and multiplied by $\varepsilon^{-1}$ to account for the surface-to-volume ratio ).

[TABLE]

Matching the orders, we obtain the following:

• $\mathcal{O}(\varepsilon^{-2})$

The first equations simply gives the independence of the leading order from $\bm{y}$ (required by homogenisation):

[TABLE]

• $\mathcal{O}(\varepsilon^{-1})$

This equation can be solved posing:

[TABLE]

Introducing the first order corrector $\bm{\chi}_{1}$, which then satisfies:

[TABLE]

• $\mathcal{O}(1)$

Matching the orders and applying Favre averaging over $\Omega_{m}$, this results in an upscaled equation for $c_{m0}$:

[TABLE]

Where $\dot{\mathcal{M}}_{i}=\int_{\partial\Omega_{i}}\bm{f}_{i}\cdot\bm{n}_{i}\text{d}S$ is the same as equation 14. Furthermore, $\bm{U}=\int_{\Omega_{m}}(\bm{u}/V_{m})\text{d}V$.

This also allows us to obtain an expression for the effective diffusivity, and to connect volume averaging and homogenisation, as $\widetilde{c_{m}}=c_{m0}$. This formulation can be seen as an alternative to spatial filtering for the mobile region, although it is possible to obtain the same correctors equations using the volume averaging method (see for example [48]).