Generalized multirate models for conjugate transfer in heterogeneous materials

We propose a novel macroscopic model for conjugate heat and mass transfer between a mobile region , where advective transport is signiﬁcant, and a set of immobile regions where diffusive transport is dominant. Applying a spatial averaging operator to the microscopic equations, we obtain a multicontinuum 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 spectral decomposition, we derive a set of equations that can be viewed as a generalization of the multirate mass transfer (MRMT) model. 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 generalized multirate transfer model (GMRT) has the advantage of providing a direct method for calculating the model coefﬁcients for immobile regions of arbitrary shapes, through the solution of appropriate microscale cell problems. An important ﬁnding is that a simple rescaling or reparametrization of the transfer rate coefﬁcient (and thus, the memory function) is not sufﬁcient to account for the ﬂow ﬁeld in the mobile region and the resulting nonuniformity of the concentration at the interfaces between mobile and immobile regions.


I. 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 a 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 Refs. [14][15][16] for an extensive review), typical flow structures and exchange phenomena arising from strong heterogeneities (see, for example, Refs. [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] cal simulations (for example, Ref. [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 (1) integrodifferential 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 multiRate mass transfer [1], which consist in modeling the transfer between mobile and immobile regions as a system of first order reactions; and (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 [13,23,24] and a unified formulation based on the multirate mass transfer has been proposed [24]. This somewhat arbitrary choice was based on the sound basis that (i) the multirate 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 multirate mass transfer can be derived from the microscopic equations exactly, and intuitively interfaced with results from homogenisation (see Refs. [25,26] for an extensive review of homogenisation theory).
However, accurate estimation for the closure parameters of the multirate mass transfer model is still a largely debated topic. Specifically, the multirate mass transfer model of Haggerty and Gorelick [1] requires a couple of parameters for each first-order reaction: (1) α HG : the apparent exchange rate coefficient and (2) β 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 and 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 [6,37] and data fitting are often employed to obtain model parameters in practice. As a result, current mathematical formulations of multirate 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 multirate 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; and (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 Sec. II, we describe the microscopic equations and the approximations we employ. In Sec. III, we present the upscaling methodology in details and in Sec. IV, we show how the model of Haggerty and Gorelick can be obtained as a zero-order approximation of our model. In Sec. IV, we also present higher-order models and we summarize the model parameters in Sec. VI. We conclude in Sec. VII with an outlook to future extensions of the current model.
Additional details on the homogenisation procedures can be found in Appendix.

A. 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 = where 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 ∂ i . We also assume that the regions i are completely included inside and that they are disconnected (i.e., they only border on m ).
In the following, we will assume that transport within inclusions i is dominated by diffusion, while on m , advection might not be negligible. Thus we can define a Peclet number: where U is a characteristic velocity of the fluid, L is a characteristic length and D is a diffusion coefficient. We will therefore assume that in the immobile regions: while no assumption is made on the Peclet number in the mobile region.

B. Microscopic governing equations
We assume that the concentration field c m (x, t ) in the mobile region is obeying the advection-diffusion equation at the microscopic scale: Furthermore, we have N i diffusion equations for the concentrations in the immobile regions: We assume here the immobile diffusion coefficients D i to be constant. This can be easily relaxed to smooth or piecewise smooth coefficients and will be subject of future studies (by decomposing into coupled subregions).
At the interfaces ∂ i , we enforce continuity of fields and fluxes: 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.

A. Spatial filtering
Standard multicontinuum models [38] can be obtained from Eqs. (3) and (4), by applying a spatial filtering operator to the governing equation for the mobile region using a REV (representative elementary volume) as support. We will assume that such REVs have a local periodic behavior or, in other words, that their geometry changes very slowly with x. Specifically, we assume that the number and geometrical configuration of the immobile regions included in a region (x) centered on x is essentially equivalent to that of a region (x + δx) for a sufficiently small δx. This procedure produces fields that are much smoother than the original ones.
It is important to notice that in our model described by Eqs. (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 macro given by the union of a number of REVs . 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 x). Therefore it is possible to interchange between and macro when computing averages without any loss of generality.
Thus we define the volume average of c m in the region (x) (but it can be trivially extended to macro ) as the top-hat filter of volume V = dV centered on x: where K 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]: where (x) is a REV centered on x. In this formulation, both c m and c m are both functions of the spatial coordinate x. However, the integral operator results in c m to be much smoother than c m and therefore we will consider as c m does not depend on space at scales smaller than . Similarly, V is also a slowly varying function of x. We can therefore write an explicit expression for c m : As commonly done for multiphase and compressible flows, we also define the Favre top-hat Kernel [39] as the volume And therefore, the Favre averaged concentration c m can be written as The spatial filtering procedure is illustrated in Fig. 2, where the resulting Favre averaged concentration is much smoother than the original one.
Generally, is itself a function of the coordinate x but, since we assume symmetry of under translation. the integral commutes with spatial derivatives and we will therefore omit its dependence of x for the sake of brevity. such that we obtain the relation: where we introduced the capacity of the mobile region β m defined as V m /V . Some authors (for example, Ref. [1]) define β m as being multiplied by a retardation factor obtained from a rescaling 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 Eq. (11), and thus to the definition of capacity.

B. Multicontinuum formulation
Assuming that the immobile regions i are fully included into (i.e., ∂ ∩ ∂ i = ∅ i = 1, . . . , N i ), applying the integral operator (6) to (3) and making use of the Green's theorem, we obtain: where we have defined the total average flux in the mobile region: and the average inter-region mass exchange rate for region i: Since in this work our focus is on the interface exchange, we do not perform an accurate upscaling of J m , and we will use a simplified expression without any loss of generality: where u eff and D m,eff are the effective velocity and the effective diffusivity in the mobile region. The capacity β 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 where V i = i dV is the volume occupied by region i . Thus we integrate (4) to obtain where β 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 x at the macro scale (for example, due to the distribution of immobile regions at the macroscale). Equation (17), substituted into (12), leads to the multicontinuum equation for the concentration field in the mobile region: In (18), we transformed the boundary conditions of the microscopic equation into source terms, one for each immobile region. However, in this formulation 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.

C. Multirate mass transfer
In order to express 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: where the function ψ i satisfies the following equation and boundary conditions: Summing Eqs. (20) and (21), and using decomposition (19) gives back Eq. (4) with the correct boundary conditions. Due to Eq. (20), ψ i satisfies the following Gauß-Green integral: being n a field normal to ∂ i and S i the surface of ∂ i . Therefore, in our formulation, the function ψ i is simply required to satisfy the nonhomogeneous time dependent boundary conditions, while c i satisfies a nonhomogeneous 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: where c i j (t ) are series coefficients that depend on time and on x at the macro scale only, while the eigenfunctions φ i j (x) carry the dependence on the spatial coordinate at the microscale and satisfy where λ i j is the eigenvalue corresponding to eigenfunction φ i j . 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 . 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 i j is only keeping track of the different initial conditions at the macroscopic scale (since the boundary conditions are homogeneous). Therefore we will take c i j out of any spatial derivative within the immobile regions, since they are assumed to be negligible.
Both Eqs. (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 D i . As a consequence, we can relate the dimensional eigenvalue λ i j with a dimensionless eigenvalue: Following this rescaling, Eqs. (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 φ ik , we obtain where A i = i φ ik φ ik dV 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: Substituting into Eq. (26), we obtain where we introduced the projection of ψ into φ ik scaled over the norm of φ ik : Therefore c i is now given by where is the correction function for the immobile region, which accounts for the nonhomogeneity of c m (x, t ) at the interface.

Computation of the exchange rate
We now compute the Favre averaged concentration in the i-immobile region: where the Favre averaged correction function is given by The terms β ik play the role of capacities (or a normalized weighting function) since Equation (35) is the so-called partition of unity (notice that β 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: It is important to notice that all the terms involved in the multirate transfer can be computed a priori by solving a cell problem, which consists in solving Eq. (20) for ψ i and the eigenvalue problem (24) for each immobile region i. However, ψ i (x, t ) is a non trivial function of c m (x, t ), 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 (x, t ) in order to make any further progress.

D. Representation of c m (x, t ) 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 (x, t ).
In order to give a representation of the spatial variability of c m (x, t ) without having to solve the microscopic unsteady governing equations, one can employ the classical two-scale expansion of homogenisation theory [26,42], and express c m as where χ n is the corrector function corresponding to the n order of the series and : represents the contraction between the corrector and the n-order gradient ∇ n . c m and ∇ n c m varies much slower than χ n in x and can be considered as constant 1 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 firstorder corrector χ 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 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 Appendix.
As a side note, we mention that homogenisation can be also employed to obtain an expression for J m,eff [43,44] 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 [45].
To introduce the information provided by the corrector equation into our problem, we can expand ψ i in a similar fashion: where the functions n are coupled with the correctors χ n at the interface and satisfy [owing the linearity of equation (20)]: These are a set of partial differential equations for tensors of rank n. Notice that in satisfies a boundary integral relation similar to Eq. (22). We can now substitute expansion (38) into where we introduced 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 (x, t ) at the interface. This formulation shows that, when we assume c m (x, t ) = c m at the interface, then θ i = 0, which means that no correction is necessary. Expansion (38) can now be substituted into the evolution equation for c ik , leading to

IV. GOVERNING EQUATIONS OF THE GENERALIZED MULTIRATE TRANSFER MODEL
We can now write down a set of equations for the generalized multirate 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 λ ik = λ ik D i /L 2 i ). The full system of equations is 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 (x, t ) 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 c m is sufficiently regular it is possible to obtain a good approximation just with the first-order corrector χ 1 [26].
A key feature of the current formulation is that accounts for the nonuniform 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 λ 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 highorder 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.

A. A note on the truncation of the multirate series: Equilibrium modes
Clearly, practical applications require the multirate series to be truncated at some value k max corresponding to β ik max , λ ik max , c ik max , and ik max n . While previous works enforced k β 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 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 max is chosen sufficiently large, according to the physics of the problem. Therefore one can write Furthermore, the corresponding eigenfunctions φ ik will be highly oscillating for k > k max in contrast with the corrector in , which we choose to be a smooth function (recall equation 20). Therefore the projection of in over φ ik will be very small for k > k max since in will not have significantly high modes. The regularity of in is also connected to the change of c m (x, t ) over ∂ i . In most of the applications (e.g., forced convection) c m (x, t ) varies regularly over the interfaces, and thus the correctors χ n vary smoothly (and slowly) over ∂ i . Therefore we can also assume that there exist a k max such that ikn = 0 fork > k max (44) and, consequently, 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 β i (and thus β m ) to account for the removed modes compared to simply rescaling β ik . In practice, one would then have new truncated immobile capacities β tc i and mobile capacities β tc m defined as 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 max in the immobile regions is completely determined by the mobile region.

B. The multirate model of Haggerty and Gorelick:
The leading order approximation The original multirate mass transfer (MRMT) model proposed by Haggerty and 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 ψ i , which results in the system: 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 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.

C. Computation of β ik and λ ik
Coefficients β ik and λ ik do not depend in any way on the interface concentration c m (x, t ) 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 I. Evaluation of λ ik and β ik for simple geometries for which there is an analytical solution of the unsteady diffusion equation (see Ref. [46] for details). Here, L represents half the length of the layer (the domain in the layer goes from −L to L), R is the radius of the sphere or cylinder and ζ k is the kth zero of the zero-order Bessel function of the first kind. Table I shows the expression of λ ik and β ik for a set of simple geometries. Clearly, our coefficients match those proposed by Haggerty and Gorelick [1], except for a factor β m in β ik , which is consistent with our formulation since we do not divide the equation for c m by β m .
In this approximation. coefficients λ ik and β ik have the same meaning as in Haggerty and Gorelick, where λ ik plays the role of exchange rate between c ik and ψ ik . As demonstrated in Table I, λ ik is a function of geometrical dimensions and material properties through λ ik = λ ik D i /L 2 i , where the dimensionless eigenvalue depends on the shape of the immobile region only. On the contrary, β ik is a dimensionless weight that depends only on the class of geometrical shapes.

V. 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: (48) Integrating the eigenvalue Eq. (24) over i , we can obtain the following relation for the eigenvalues: Expanding the first term on the right-hand side of Eq. (48) leads to where we employed Eq. (22) on the right-hand side. Then, substituting expansion (38) results into 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.

A. Generalized multirate transfer equations
The complete set of equations (42) can be therefore rewritten without mixed time-space derivatives as where we defined the equilibrium concentration for term k of region i as 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 ∂c ik /∂t = 0.

B. 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 multirate equation for the mobile region. The governing equations are given by For the special case in which the material microstructure does not vary in space and the flow field is macroscopically homogeneous (i.e., u eff = const), ik1 does not depends on the spatial coordinates and we can define a drift velocity: and thus a new effective velocity Therefore the equation for the mobile region simply reduces to a standard advection diffusion equation, with an additional multirate reactive term:

C. Physical considerations on ik1
Equation (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 nonequilibrium at the microscale (i.e., within ) and can be attributed (at least asymptotically) to 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 nonequilibrium.
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 characterized by having Now, we consider Eq. (37) at the first order: 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 χ 1,x : Therefore, if all components of ∇ c m are positive, all the components of χ 1 are also positive. This positivity is transferred to i1 through Eq. (39) due to the properties of elliptic operators. While ik1 is not necessarily positive, the projection on the first eigenfunction i11 is positive.
It is easy to demonstrate that this positivity property holds also when ∇ c m < 0.
Therefore, as illustrated in Fig. 3, a in a system with ∇ c m > 0, the equilibrium concentration in the immobile regions will be larger than c m due to the higher value of c m at the interface. On the contrary, in the case ∇ 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.

D. 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): Now, we can decompose tensor ik2 into hydrostatic and deviatoric components: where tr( ik2 ) is the trace of ik2 and I is the identity tensor. We now introduce the diffusion coefficient arising from the conjugate transfer D ct : 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 D 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: Therefore the equation for the mobile concentration simplifies to

VI. SUMMARY OF MODEL PARAMETERS
Clearly, the multirate 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; [λ ik ]: these are simply the eigenvalues corresponding to the homogeneous eigenproblem in the immobile region, and β ik : these weights can be calculated similarly to λ ik , from the solution of the eigenproblem. Once the eigenfunctions are known, β ik is given by These are the same parameters of standard multirate 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 χ in obtained from homogenisation theory.
ikn . Projection of the function in on the eigenfunction φ ik scaled with the norm of φ ik . Clearly, the number of these parameters equals the number of terms in the multirate expansion but one can exploit some knowledge of the microstructure to simplify their expression.
Other quantities we introduced, like D ct , c eq ik , or u 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 λ ik , β ik , and ikn become irrelevant, it is still crucial to remember that all the physics of non equilibrium is contained in ikn .

VII. CONCLUSIONS
In this paper, we propose a novel approach to derive the multirate mass transfer model that is different from that of the memory function or that of Haggerty and 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 and Gorelick, it also contains their multirate model as a special case and allows extension to nonequilibrium 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 multirate 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 modeling real systems. Some examples are exchange between immobile regions; multiple mobile regions with different mobility (e.g., fractures); chemical reactions at interfaces; multiphase flow, heat, and mass transfer; and 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.

ACKNOWLEDGMENTS
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 [47]."

APPENDIX: 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 Eq. (12) can be obtained using homogenisation and how to calculate the first-order corrector χ 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, Ref. [44].
The process starts defining an expansion parameter and a microscopic scale y such that where R is a characteristic length of the immobile regions and L is a characteristic length at the macroscale. The field c m (x, y, t ) is then expanded in asymptotic series of ε: c m (x, y, t ) = ∞ n=0 ε n c mn (x, y, t ).
(A3) Spatial differential operators are expanded to account for the microscopic scale: Here, y represents the variation across the REV , while x is a coordinate on the macroscopic volume macro . This splitting is also known as two-scale asymptotics. The boundary conditions are of the second type, in agreement with Eq. (5): where the order ε is taken due to the scaling of 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 ε −1 to account for the surface-to-volume ratio).
Matching the orders, we obtain the following. O(ε −2 ). The first equations simply gives the independence of the leading order from y (required by homogenisation): c m0 = c m0 (x, t ).
whereṀ i = ∂ i f i · n i dS is the same as Eq. (14). Furthermore, U = m (u/V m )dV . This also allows us to obtain an expression for the effective diffusivity, and to connect volume averaging and homogenisation, as 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, Ref. [48]).