Tristan Babey
email: tristan.babey@univ-rennes1.fr

Jean-Raynald De Dreuzy

Céline Casenave

Multi-Rate Mass Transfer (MRMT) models for general diffusive porosity structures

Keywords: Porous media, Solute transport, Mobile-immobile models, Multi-Rate Mass Transfer

. We develop effective shapefree numerical methods to identify its few dominant rates, that comply with any distribution of rates and porosities. We show that differences in terms of macrodispersion are not larger than 50% for approximate MRMT models with only one rate (double porosity models), and drop down to less than 0.1% for five rates MRMT models. Low-dimensional MRMT models BABEY ET AL.: MRMT MODELS FOR GENERAL DIFFUSIVE POROSITY STRUCTURES 2 accurately approach transport in structured diffusive porosities at intermediate and long times and only miss early responses.

Introduction

Transport in complex geological environments results in part from the interactions between fast advective-dominated transport in a localized "mobile porosity" and slow diffusivedominated transport in extensive "immobile porosities". It is the case of the fracture-matrix systems [START_REF] Neretnieks | Diffusion in the Rock Matrix: An Important Factor in Radionucleside Retardation?[END_REF][START_REF] Tang | Contaminant transport in fractured porous-mediaanalytical solution for a single fracture[END_REF] and of the highly heterogeneous porous media [START_REF] Fernandez-Garcia | Upscaling transport with mass transfer models: Mean behavior and propagation of uncertainty[END_REF][START_REF] Golfier | Comparison of theory and experiment for solute transport in highly heterogeneous porous medium[END_REF][START_REF] Gotovac | Flow and travel time statistics in highly heterogeneous porous media[END_REF][START_REF] Willmann | Transport upscaling in heterogeneous aquifers: What physical parameters control memory functions?[END_REF]. When diffusive times in the immobile zones become much larger than the characteristic advective time in the mobile zone, transport becomes anomalous with non-Gaussian concentration plumes, more extensive spreading and mixing, slow transit times, and broad ranges of solute retardation times [START_REF] Berkowitz | Modeling non-Fickian transport in geological formations as a continuous time random walk[END_REF][START_REF] Dentz | Time behavior of solute transport in heterogeneous media: transition from anomalous to normal transport[END_REF]. Such transport mechanisms and exchanges are at the root of numerous anomalous transport modeling frameworks [START_REF] Benson | A simple and efficient random walk solution of multirate mobile/immobile mass transport equations[END_REF][START_REF] Benson | Application of a fractional advection-dispersion equation[END_REF][START_REF] Berkowitz | Theory of anomalous chemical transport in random fracture networks[END_REF][START_REF] Carrera | On matrix diffusion: formulations, solution methods and qualitative effects[END_REF][START_REF] Cushman | Fractional advection-dispersion equation: A classical mass balance with convolution-Fickian flux[END_REF][START_REF] Haggerty | Multiple-Rate Mass Transfer for Modeling Diffusion and Surface Reactions in Media with Pore-Scale Heterogeneity[END_REF] and can be highly effective in the interpretative and predictive phases of laboratory and field experiments [START_REF] Benson | Fractional dispersion, Lévy motion, and the MADE tracer test[END_REF][START_REF] Berkowitz | Anomalous transport in laboratory-scale, heterogeneous porous media[END_REF][START_REF] Gouze | Non-Fickian dispersion in porous media explained by heterogeneous microscale matrix diffusion[END_REF][START_REF] Haggerty | Tracer tests in a fractured dolomite 2. Analysis of mass transfer in single-well injection-withdrawal tests[END_REF][START_REF] Haggerty | What controls the apparent timescale of solute mass transfer in aquifers and soils? A comparison of experimental results[END_REF]. Anomalous transport ultimately stems from some extended distribution whether it is a waiting time distribution as in Continuous Time Random Walk (CTRW) or a rate-porosity distribution as in Mutli-Rate Mass Transfert (MRMT) [START_REF] Dentz | Transport behavior of a passive solute in continuous time random walks and multirate mass transfer[END_REF][START_REF] Neuman | Perspective on theories of non-Fickian transport in heterogeneous media[END_REF][START_REF] Silva | A general real-time formulation for multi-rate mass transfer problems[END_REF]. For MRMT models, while these distributions can take very different shapes [START_REF] Haggerty | On the late-time behavior of tracer test breakthrough curves[END_REF], only some power-law distributions are effectively related to diffusive processes in 1D, 2D or 3D inclusions [START_REF] Carrera | On matrix diffusion: formulations, solution methods and qualitative effects[END_REF][START_REF] Haggerty | Multiple-Rate Mass Transfer for Modeling Diffusion and Surface Reactions in Media with Pore-Scale Heterogeneity[END_REF] or to anomalous diffusive processes in fractal-like structures [START_REF] Haggerty | Matrix diffusion -heavy-tailed residence time distributions and their influence on radionuclide retention[END_REF]. Diffusive structures may however be topologically more complex like for example for fracture dead ends [START_REF] Flekkøy | Modeling hydrofracture[END_REF][START_REF] Sornette | Fault Growth in Brittle-Ductile Experiments and the Mechanics of Continental Collisions[END_REF], fracture-matrix interactions [START_REF] Jardine | Quantifying diffusive mass transfer in fractured shale bedrock[END_REF][START_REF] Karimi-Fard | Generation of coarse-scale continuum flow models from detailed fracture characterizations[END_REF][START_REF] Sudicky | Contaminant transport in fractured porous-media -analytical solutions for a system of parallel fractures[END_REF][START_REF] Tang | Contaminant transport in fractured porous-mediaanalytical solution for a single fracture[END_REF][START_REF] Tsang | Study of alternative tracer tests in characterizing transport in fractured rocks[END_REF], or dissolution patterns in porous media [START_REF] Golfier | On the ability of a Darcy-scale model to capture wormhole formation during the dissolution of a porous medium[END_REF][START_REF] Luquot | Experimental Characterization of Porosity Structure and Transport Property Changes in Limestone Undergoing Different Dissolution Regimes[END_REF] (Figure 1).

In this article, we show that the MRMT framework is general to all diffusive architectures that can be modeled as a finite number of interconnected continua (Figure 1). The notion of continuum comes from the double porosity and Multiple INteracting Continua (MINC) concepts introduced initially for fracture-matrix systems [START_REF] Pruess | A practical method for modeling fluid and heat-flow in fractured porous-media[END_REF][START_REF] Warren | The Behavior of Naturally Fractured Reservoirs[END_REF]. The double porosity model is the classical diffusive interaction of advective-diffusive processes in a mobile zone with a single immobile zone like in double-porosity models [START_REF] Warren | The Behavior of Naturally Fractured Reservoirs[END_REF]. The Multiple INteracting Continua (MINC) framework models matrix diffusion as diffusive-like exchanges within a succession of "continua", identified to the elementary cells issued by a finite-difference discretization of the diffusion process in the matrix (Figure 2a) [START_REF] Pruess | Brief Guide to the MINC -Method for Modeling Flow and Transport in Fractured MediaRep[END_REF][START_REF] Pruess | A practical method for modeling fluid and heat-flow in fractured porous-media[END_REF]. The denomination of multiple continua is a direct generalization of the double porosity concept of [START_REF] Warren | The Behavior of Naturally Fractured Reservoirs[END_REF]. We propose to further generalize the notion of interacting continua to any immobile zones structure where diffusive-like exchanges intervene between any connected zones or continua. Because of the potential importance of structure on diffusion, we denote these models as Structured INteracting Continua (SINC). SINC models include a wide range of structures going from elementary branching and loops (Figure 2b andc) to more involved dissolution patterns (Figure 2d). They would typically be derived from the coarse discretization of diffusion processes in dead-end porosity structures [START_REF] Gouze | Non-Fickian dispersion in porous media explained by heterogeneous microscale matrix diffusion[END_REF][START_REF] Noetinger | Up-scaling of double porosity fractured media using continuous-time random walks methods[END_REF]. We define SINC models in section 2, with their exact relation to the MRMT and MINC models. We show in section 3 that any SINC model is equivalent in terms of transport to a unique MRMT model of the same dimension, i.e. with the same number of immobile zones. We develop efficient numerical methods in section 4 to identify lower-dimension but highly accurate approximate MRMT models. The SINC model is made up of a continuous 1-D mobile zone in interaction with a finite number of interconnected immobile zones (Figure 1). In the continuous mobile zone, solutes are transported by advection, dispersion and diffusion. Between the mobile zone and the immobile zones as well as between the immobile zones, concentration exchanges are diffusive-like, i.e. directly proportional to the difference of concentrations. This model can be generically expressed as:

  U R L AU t U m     (1)
where

T N im im m t x c t x c t x c U )] , ( ) , ( ) , ( [ 1  
is the vector of dimension N+1 made up of the concentrations in the mobile zone

  t x c m
, and in the N immobile zones

  t x c i im , with N i ,..., 1  . A is the ) 1 , 1 (   N N
interaction matrix characterizing the diffusive-like concentration exchanges between the immobile zones and with the mobile zone. R m is the restriction matrix to the mobile zone:

      1 1 ,    j i j i R m   . ( 2
)
L is the transport operator in the mobile zone:

  2 2 x c d x c q c L m m m m m         (3)
with q, m  and m d the Darcy flow, porosity and dispersion coefficient in the mobile zone.

The physical properties of the mobile and immobile domains are homogeneous along the mobile domain. The interaction matrix A is equal to the matrix M deriving from the diffusivelike mass exchanges between the different zones, corrected by their porosities:

M A 1     (4)
with  the diagonal matrix whose diagonal elements are the porosities associated to the mobile zone  m and to the N immobile zones

i im  with N i ,..., 1  :                    1 , 1 , for , , 1 , 1 1 j i j i j i i im m    . ( 5
)
The matrix M expresses rates of mass exchange. The dimension of its elements is therefore the inverse of a time. As exchanges are diffusive-like, M is a M-matrix, i.e. M is symmetric, its diagonal elements are positive, its off-diagonal elements are negative or equal to zero, and the sum of its elements along each of its rows is equal to zero. We underline that it is the matrix M that is symmetric and generally not the interaction matrix A that also integrates the differences in porosities. The interaction matrix A registers the connectivity of the different zones through the position of its non-zero off-diagonal elements, the strength of the interactions is determined by porosity ratios and exchange rates. Figure 3a shows the example of the interaction matrix for the asymmetric Y structure (cf Figure 2b). The branching architecture leads to a compact interaction matrix, which values are all on the three principal diagonals except at the branching node. The interaction matrix is scaled by the inverse of the mean diffusion time in the immobile structure . We define  as the quadratic mean distance of the immobile zones to the mobile zone divided by the diffusion coefficient between two zones.

The SINC framework generalizes Multiple INteracting Continua (MINC) [START_REF] Pruess | A practical method for modeling fluid and heat-flow in fractured porous-media[END_REF]. MINC models are obtained by the finite-difference discretization of diffusion in a 1-D homogeneous medium. In such a case, the interaction matrix A is simply tri-diagonal. The term "continuum" comes from the continuity along the mobile zone like in the dual-porosity concept [START_REF] Warren | The Behavior of Naturally Fractured Reservoirs[END_REF]. It is consistent with our concept of a continuous formalism of transport along the mobile zone and a finite number of interacting continua in the immobile direction.

The SINC framework also generalizes Multiple-Rate Mass Transfer models with a finite number of rates [START_REF] Carrera | On matrix diffusion: formulations, solution methods and qualitative effects[END_REF][START_REF] Haggerty | Multiple-Rate Mass Transfer for Modeling Diffusion and Surface Reactions in Media with Pore-Scale Heterogeneity[END_REF]. In MRMT models as defined by [START_REF] Haggerty | Multiple-Rate Mass Transfer for Modeling Diffusion and Surface Reactions in Media with Pore-Scale Heterogeneity[END_REF], the immobile domain consists in a distribution of subdomains exchanging exclusively with the mobile domain. The star-like connectivity structure leads to an arrow-type broad-width interaction matrix (Figure 3b). Each sub-domain is characterized by its rate of transfer i  and its porosity i im  , with [START_REF] Haggerty | Multiple-Rate Mass Transfer for Modeling Diffusion and Surface Reactions in Media with Pore-Scale Heterogeneity[END_REF]. When defined with a finite number of rates, MRMT models can be recovered by fixing

1   i i   [
                                              ) 1 , 1 ( ) , ( for , , 1 , 1 , , 1 for , 1 1 , and 1 , 1 for 0 , 1 , 1 1 j i j i j i j i M i i M i i M i M j i j i j i M i im m i j j i i im      . ( 6
)
Because they are defined in algebraic terms, SINC models can only represent MRMT models with finite number of rates. They are rigorously not equivalent to 1D, 2D and 3D diffusion models that themselves correspond to an infinite number of exchange rates [START_REF] Haggerty | Multiple-Rate Mass Transfer for Modeling Diffusion and Surface Reactions in Media with Pore-Scale Heterogeneity[END_REF]. However they offer very accurate approximations when taking only a finite number of rates, especially as the porosities i im  decrease as a power law of i. SINC models generalize MRMT models defined with a finite number of rates but remain different to MRMT models defined by an infinite serie of rates or by a continuous rate-porosity function.

Our objective in this article is not to produce any kind of generalization of MRMT models but rather to determine how complex diffusive porosity structures can be approached by simpler reduced models. 1 and the multiplicative factor  (=5.015) is equal to the ratio of the distance between two consecutive immobile zones to the radius of gyration of the immobile domain to the mobile zone.  is the ratio of the total immobile porosity to the mobile porosity.

(c)

(d)
Following up the work of [START_REF] Haggerty | Multiple-Rate Mass Transfer for Modeling Diffusion and Surface Reactions in Media with Pore-Scale Heterogeneity[END_REF], we define the equivalence relation to a MRMT model as the identity of the mobile concentration at any time. To prove the equivalence of the SINC model to a MRMT model, we first decompose the interaction matrix A in two parts, consisting in the exchanges between the immobile zones only for the first part and in the exchanges between the mobile zone and the immobile zones for the second one:

                                1 , 1 21 1 , 1 12 11 0 0 0 0 N N A A A A A A A        (7) with ) 1 2 , 1 2 (    N N A A    . As shown in Appendix A, the matrix A  can be diagonalized 1    R R A    (8)
with  the diagonal matrix made up of the eigenvalues of A  , which are all real and negative, and R  the matrix of the eigenvectors of A  , defined each up to a multiplicative constant. To fix the eigenvalues decomposition, as  is defined up to a permutation of its diagonal elements, we sort the eigenvalues according to their absolute value by increasing order

1 , 1 ,      i i i i
. In the new coordinate system defined by the eigenvectors, there are no more exchanges between the immobile zones. All exchanges are made directly between the immobile zones and the mobile zone. To characterize the precise nature and extent of these exchanges, we propagate the change of the coordinate system to the exchanges with the mobile zone

1 1 , 1 21 1 , 1 12 11 1 , 1 21 1 , 1 12 11 0 0                                   R B B B B A R A A A A A N N N N     (9) with the ) 1 , 1 (   N N matrix R and its inverse 1  R defined as                0 0 0 0 1 R R    and                  0 0 0 0 1 1 1 R R    . ( 10
)
In this transformation, the mobile zone remains unchanged consistently with our objective to reorganize only the immobile zones without interfering with the concentration in the mobile zone. The full transformation defined by R is applied to the matrix A following its decomposition in equation ( 7)

1 1 , 1 21 1 , 1 12 11 1 0 0 0 0                                    R B B B B A R R R A N N       (11)
which can be finally expressed by a simple factorization as

1   RBR A with                   1 , 1 21 1 , 1 12 11 N N B B B B A B   . ( 12
)
Finally, as the restriction matrix to the mobile zone m R (equation ( 2)) commutes with 1  R , we can write the full model as

  U R R L U BR t U R m 1 1 1        . ( 13
)
To be representative of a MRMT model, B should be an arrow matrix, the sum of its elements over each of its rows should be zero, and all its non-diagonal elements must be either positive or equal to zero. In Appendix B, we show that this can be obtained by adjusting the norm of the eigenvectors in R  . The characteristics of this MRMT model are then determined by a simple identification of equation ( 6) with equation ( 12):

               i i m i i m i im i i B B B 1 , 1 1 , 1 1 , 1     . ( 14
)
While this algebraic identification method can already be widely used to determine equivalent MRMT models, it faces some numerical limitations. The diagonalization process becomes challenging when A ˆ becomes large, limiting the range of the identification to not too complex architectures and/or coarse discretizations of the immobile domain. The immobile porosity structure may be composed of a large number of cells, while a much smaller number of rates may be necessary to get highly accurate equivalent MRMT models. To address these limitations, we propose an alternative approximate numerical identification method.

Approximate numerical identification method of the MRMT model equivalent to a SINC model

We first develop the numerical methodology and secondly apply it to the four examples of properties, the same overall immobile to mobile porosity ratio , and the same mean quadratic diffusion time  as defined in section 2 (Table 1). The dispersivity in the mobile zone

  m m q d  / /
divided by the effective dispersivity due to the exchanges with the immobile zone defined as   m q   / is taken equal to 5 10 . 5  , i.e. much smaller than 1, so that dispersive effects come predominantly from the mobile/immobile exchanges. For the same reason,  is taken 1: Parameters used for the simulation of section 4.2 with the characteristic diffusion time  and the consecutive distance covered by advection in the mobile zone m q   / as temporal and spatial dimensional parameters.  is the immobile to mobile porosity ratio.

much larger than 1 ( 100   ). Parameter Value  100     2 / / m m q d   5 10 . 5    m q    / / 0 3.10 -3 0 m ax /  x 100 0 /  dx 0.2  / dt 4 10 . 5  Table
    2 / / m m q d   is the dimensionless dispersion in the mobile zone.   m q    / / 0 is the dimensionless standard deviation of the initial Gaussian concentration profile. m ax
x is the extension of the simulation domain in the direction of the mobile zone, dx is the spatial step along the mobile zone and dt is the time step.

Methodology

To numerically approximate the MRMT model equivalent to a SINC model, we consider the particular case of the discharge of the immobile zones to the mobile zone. The mass discharged to the mobile zone from an initially homogeneous immobile concentrations 0 c is fitted by a combination of exponential functions typical of the MRMT model. The residual mass per unitary volume m(t) in the immobile domains is fitted by its MRMT counterpart

) (t 
given by:

    N i i i im t c t 1 0 ) exp( ) (    ( 15
)
where

i  and i im  , N i  1 
are the N rates and immobile porosities defined in equation ( 6).

In Appendix C, we develop an optimization method to identify the  and im  series. We further illustrate its application to the cases

1  N and 2  N in Appendix D.
The advantage of this method over the previous algebraic method is to be flexible in terms of the number of rates N to identify. It also prioritizes the identification of the rates and immobile porosities having a significant impact on transport. In the following we determine the number of rates that should be identified to model accurately macrodispersive processes.

Simulations and results

We analyze the influence of the number of rates N on the reproduction of macrodispersion for the four structures displayed in Figure 2 for N ranging from 1 to 5. Identification is performed according to the methodology presented before and with a logarithmic sampling of times starting at the time for which the relative mass discharge to the mobile zone is lower than 10 -3 and ending at the time for which the relative residual mass itself is lower than 10 -4 .

Simulations of transport for the SINC and approximate MRMT models are further performed with identical Gaussian concentration profiles in the mobile and immobile domains

              2 0 2 0 0 max 2 ) ( exp 2 ) 0 , ( ) 0 , (    x x c t x c t x c im m (16)
where 0

x , 0  and m ax c are the mean, standard deviation and maximum concentration of the Gaussian profile. 0  is taken small enough so that the initial plume size has minor effects on the overall dispersion (Table 1).

16

We compute for all models the effective dispersion coefficient D as

dt d D x 2 2 1   (17)
with  x the plume spreading

2 0 1 0 2 2 ) / ( / m m m m x    . ( 18
)
The spatial moments of concentration k m are given by

        0 1 1 ) , , ( ) , ( ) ( x N i k k dx i t x U i i x t m . ( 19
)
We assess the quality of the MRMT model for modeling dispersion through the quadratic mean of the relative difference in effective dispersion of the SINC and MRMT models:

           t t MRMT SINC MRMT SINC MRMT SINC dt dt t D t D t D t D D D diff 2 2 / )) ( ) ( ( ) ( ) ( ) , ( . (20)
This criterion is sensitive to differences in dispersion on a broad time range from the initial time to the time at which dispersion becomes constant.

Classical numerical methods were used to solve both the exchanges within the immobile zones and the advective-dispersive transport in the mobile zone [START_REF] De Dreuzy | Influence of porosity structures on mixing-induced reactivity at chemical equilibrium in mobile/immobile Multi-Rate Mass Transfer (MRMT) and Multiple INteracting Continua (MINC) models[END_REF]. They were validated on a set of immobile structures rigorously equivalent to the MINC model (Figure 4). In those cases, ) , (

MINC SINC D D diff
was of the order of 10 -11 and much smaller than any differences recorded in the analysis performed hereafter. For the structures displayed on Figure 2, the effective dispersion obtained with the MRMT models converges quickly to the reference dispersion of the corresponding SINC model (Figure 5). The equivalent MRMT model with only one immobile zone (N=1), equivalent to the double porosity model [START_REF] Warren | The Behavior of Naturally Fractured Reservoirs[END_REF], already gives the right order of magnitude of dispersion. With N=2, the error of the MRMT model is close to10%. With N=4, it is close to 1% and with N=5, it is close to 0.1%. A very limited number of rates is thus sufficient to represent even the complex diffusive structure displayed by Figure 2d. This fundamentally comes from the homogenization nature of diffusion that systematically removes the extremes of the concentration distributions as previously noted in numerous studies [START_REF] Haggerty | Multiple-Rate Mass Transfer for Modeling Diffusion and Surface Reactions in Media with Pore-Scale Heterogeneity[END_REF][START_REF] Villermaux | Chemical-engineering approach to dynamic modeling of linear chromatography -a flexible method for representing complex phenomena from simple concepts[END_REF].

In addition, the equivalent MRMT model with only one rate reproduces well the tailing of the breakthrough curve (Figure 6). As expected, introducing higher rates progressively improves the accuracy of the MRMT model at earlier times. The double peak observed for N=1 is a classical feature of double porosity models where advection is much faster than diffusion in the immobile porosity [START_REF] Michalak | Macroscopic behavior and random-walk particle tracking of kinetically sorbing solutes[END_REF]. It vanishes for higher-order MRMT models (N=2 to 5), as higher rates enhance short-term mobile-immobile exchanges and remove early breakthroughs.

The quality of the MRMT model with only very few rates fundamentally comes from the dominating role of the smaller rates (i.e. larger transfer times). In fact the whole rate series as determined by the algebraic diagonalization method shows that the lowest rate dominates in every case by accounting for 70% to 85% of the total immobile porosity (Figure 7). The five lowest rates represent at least 95% of the total immobile porosity for all the studied structures.

The evolution of the porosities i im  with the rates i  is monotonic only for the MINC model and becomes much more variable for more complex structures, highlighting the need of identification methods that do not assume any a priori repartition of the immobile porosity among the rates. We finally note that our results pertain to complex diffusive structures observed on a given range of scales. If not close to the mobile zone, finer details would be fast homogenized by diffusion and are unlikely to modify the identified rates. If close to the mobile zone, finer details should be treated independently in the same way as for the largerscale structures and the MRMT models obtained at different scales should be eventually superposed. A simple sensitivity study on representative diffusive structures shows that very few rates are needed for accurately modeling the solute transport in a 1D advection-dominated mobile zone exchanging with an immobile porosity structure. Double porosity models (MRMT with a single rate) already give the right order of magnitude of macrodispersion. Differences in macrodispersion drop down to around 10% for two rates, to 1% for four rates, and to less than 0.1% for five rates. Simplified models based on only five rates approach accurately the behavior of the system at intermediate to large times and only miss the very early responses.

While only few rates are necessary, their distribution and associated porosities are highly variable, the complexity of the structure being transferred to the identified rates and associated porosities. We thus conclude that MRMT models can be very efficient for modeling diffusion-like transport in a broad range of porosity structures with only very few rates. Even though numerical simulations have been done in 1D mobile domains, results are likely generalize to 3D. Additional simulations should also be performed to investigate the behavior of mixing and chemical reactivity both between different SINC structures and between the SINC structures and their simplified MRMT counterparts.

Appendix A: Diagonalization of A



We show that the eigenvalues of the matrix A  are real and negative and that they correspond to the opposites of the rates of the equivalent MRMT model ( i   ). As displayed by equations ( 4) and ( 7),

M A    1     where ) 1 2 , 1 2 (      N N    and ) 1 2 , 1 2 (    N N M M    .
1    is diagonal and its diagonal elements are all positive.

Thus 1    is positive definite, i.e. for every non-zero and real column vector x,

0 ) , ( ) ( 2 1       i T i i i x x x   . M  is symmetrical, real, diagonally dominant    i j j j i M i i M , ) , ( ) , ( 

Appendix B: Construction of the matrix B

The norm of each eigenvector in R  is defined up to a constant. A first straightforward step is to adjust these norms so that the sum along the 2..N+1 rows of B is equal to zero. It is achieved by taking

i i B     1 , 1
, overly written in matrix form

                        1 1 1 , 1 1 , 2   N B B . ( B1
)
Given this choice and the properties of the matrices A and R  , we demonstrate that the sum of the elements of the first line of B is also zero. We express the relation between the first columns of A and B from equation ( 12):

                        1 , 1 1 , 2 1 1 , 1 1 , 2 N N A A R B B    . (B2)
As the sum of the elements of each line of A is zero

                       1 1 1 , 1 1 , 2    A A A N (B3) equation (B2) rewrites                                       1 1 1 1 1 1 1 , 1 1 , 2       R A R B B N . ( B4
)
By substituting equation (B4) into equation (B1), we deduce that the eigenvectors comply with

                     1 1 1 1    R . (B5) As U R 1 
corresponds to the concentrations in the equivalent model (equation ( 13)), equation (B5) implies that a homogeneous immobile concentration profile in SINC remains unchanged in the equivalent model. We finally express the sum of the 2...N+1 elements of the first row of B and use the result of equation(B5):

                                                   1 2 1 , 1 , 1 1 , 1 2 , 1 1 , 1 2 , 1 1 , 1 2 , 1 1 2 , 1 1 1 1 1 1 1 N j j N N N N j j A A A A R A A B B B        . ( B6
)
An additional condition for B to be representative of a MRMT model is

0 , 1  j B for 1 2...N+ j =
. In the following we show that the adjustment of the norms of the eigenvectors is sufficient to ensure this condition. Equations ( 12) and (B2) give:

              R R B B R A A R A A M R M M R M M R A A B B N im im T N m N im im N m N N im im m N m N m N N                                                                 1 1 , 1 1 , 2 1 1 , 1 1 , 2 1 , 1 1 , 2 1 1 , 1 1 , 2 1 , 1 2 , 1 1 , 1 2 , 1 1 , 1 2 , 1 1 1 1 symmetric is because 1 1 . (B7) We now show that the matrix R R T     with ) 1 2 , 1 2 (      N N   
is diagonal with only positive diagonal elements.

We note

M A    1     with ) 1 2 , 1 2 (    N N M M   
. M  is symmetric, its diagonal elements are positive, its non-zero off-diagonal elements are negative, but the sum of its elements over each of its rows is not equal to zero. As the i im  are positive, we can consider the root of the matrix   :

2 / 1 2 / 1 2 / 1 2 / 1 C ) 1 (                          M A . ( B8
) C 
is similar (in the mathematical sense) to A  , so it is diagonalizable and has the same eigenvalues as

A  . Moreover C 
is symmetric, so it can be diagonalized by an orthogonal matrix S:

1    S S C  (B9)
As a consequence, equation (B8) rewrites:

1 2 / 1 2 / 1 2 / 1 1 2 / 1 ) (              S S S S A      . (B10)
As the norms of the eigenvectors in R  are adjusted so the equation (B1) is verified, there exists a unique orthogonal matrix S such that

R S   2 / 1   . ( B11
)
As S is othogonal, S T S is diagonal with only positive diagonal elements and writes

R R R R R R S S T T T T                    2 / 1 2 / 1 2 / 1 2 / 1 ) ( . (B12) The matrix R R T   
 , which is present in equation (B7), is thus diagonal with only positive diagonal elements. Consequently, as the immobile zone and when the immobile zones are constantly discharging to a quickly flushed mobile zone where concentration is assumed to remain negligible ( [START_REF] Haggerty | Multiple-Rate Mass Transfer for Modeling Diffusion and Surface Reactions in Media with Pore-Scale Heterogeneity[END_REF], Appendix B). In the Laplace domain, exponential functions become simple rational functions and

) ( ) ( ~p p m   is expressed as     N i i i i im p c p m 1 0 1 1 ) ( ~   (C1)
where p is the Laplace variable and ) ( ~p m (respectively

) ( ~p  ) is Laplace transform of m(t) (respectively ) ( p 
). We multiply equation (C1) by the polynomial P(p) of degree

N N N N i i p a p a p a p p P          ... 1 ) 1 ( ) ( 2 2 1 1  (C2) and obtain                      N i N i j j j i i im N N p c p m p a p a p a 1 1 0 2 2 1 ) 1 ( ) ( ) ... 1 (    . ( C3
)
If we now consider the polynomial

) ( p Q i of degree N-1 1 1 , 2 2 , 1 , 1 ... 1 ) 1 ( ) (             N N i i i N i j j j i p a p a p a p p Q  (C4)
and substitute equation (C4) into equation (C3), we obtain: (C5)

) ... 1 ( ) ( ) ( ) ... ( 1 1 1 , 2 2 , 1 , 0 2 2 1              N i N N i i i i i im
The interest of equation ( C5) is to be linear in the a i (polynomial coefficients of P(p)) and in

i im  / i  with N i  1 
. We isolate these quantities from the Laplace parameter-dependant elements to obtain the linear system

                                                                                     . ( C7
                              N N N N im im im G 2 2 1 1 2 1         (C8) with                N N N N N N N a c a c a c a c c c G     1 , 0 1 1 , 1 0 1 , 0 1 1 , 1 0 0 1 0      (C9)
where the values of of p by using a least-square method, minimizing the mismatch objective function J [START_REF] Garnier | Direct Identification of Continuous-time Models from Sampled Data: Issues, Basic Solutions and Relevance, in Identification of Continuous-time Models from Sampled Data[END_REF][START_REF] Ljung | System identification : theory for the user[END_REF]:

  2 1 ) ( ) (     K k k T k p p y J   . ( C10
)
The sampling should be extensive enough to contain all the information necessary to identify the different rates. If The minimum  ~ of J is explicitly given by:  ~can also be obtained in the temporal domain. We first divide equation (C5) by N p (which is equivalent to integrate N times over t in the temporal domain) because of the better numerical stability of integration compared to derivation

) ( ) ( ) ( ) ( ~1 1 1 k K k k k T K k k p y p p p                 (C11)
) 1 ... 1 1 1 ( ) ( ) ... 1 1 1 ( 1 1 , 2 2 , 1 1 , 0 2 2 1 1 p a p a p a p c p m a p a p a p N i N i N i N i N i i im N N N N                   . (C12)
For convenience, we note

1 2 0 0 0 0 ) ( 1 2 1 ) ( ) ( du du du u m du u m t u u u n n n n          . ( C13
)
The inverse Laplace transform of equation (C12) gives

) ... )! 3 ( )! 2 ( )! 1 ( ( ) ( ... ) ( ) ( ) ( 1 1 , 3 2 , 2 1 , 1 0 ) 2 ( 2 ) 1 ( 1 ) (                        N i N i N i N i N i i im N N N N a N t a N
N N i N i i i im N i i i i im N i i i im N N N N N                                                                                        ) ( 1 1 , 0 1 1 , 0 1 0 2 1 2 1 ) 2 ( ) 1 ( ) ( ) ( , , 1 )! 2 ( )! 1 ( ) ( ) ( ) ( ) (             . ( C15
)
 may also be obtained with a similar least-square method by considering a discretization of time k t (k=1...K) and by minimizing the objective function

  2 1 ) ( ) (     K k k T k t t y J   . (C16)
The minimum  ~ is given by Dividing equation (D3) by p and then using the inverse Laplace transform, we obtain:

) ( ) ( ) ( ) ( ~1 1 1 k K k k k T K k k t y t t t                 . ( C17
1 0 1 1 0 ) ( ) (    c t m du u m im t    (D4)
Equations ( D3) and (D4) are both linear in quantities depending on the unknown parameters 1  and 1 im  to be identified. Equation (D4) can be written under the form:

) ( ) ( t y t T    (D5) with:         1 ) ( ) ( t m t  ,        1 1 0 1 / / 1     im c ,    t du u m t y 0 ) ( ) ( . ( D6
)
The vector  ~ which minimize the quantity:

  2 1 ) ( ) (     K k k T k t t y J   (D7) with a discretization k t , K k ,..., 1  of t, is given by: ) ( ) ( ) ( ) ( ~1 1 1 k K k k k T K k k t y t t t                 (D8) where          1 ) ( ) ( ) ( ) ( ) ( 2 k k k k T k t m t m t m t t   ,                   k k t t k k k du u m du u m t m t y t 0 0 ) ( ) ( ) ( ) ( ) (  (D9) so                                        K k t K k t k K k k K k k K k k k k du u m du u m t m K t m t m t m 1 0 1 0 1 1 1 1 2 ) ( ) ( ) ( ) ( ) ( ) (  . (D10)
We then get 1  and 1 im  :

1 1 1    , 0 1 2 1 ~c im     . (D11) Case 2  N In Laplace domain, equation (D1) rewrites: 2 2 0 2 1 1 0 1 1 1 1 1 ) ( ~      p c p c p m im im     (D12) which gives when multiplied by          . (D13)    (D16) with:                     1 ) ( ) ( ) ( 0 t t m du u m t t  ,                       2 1 0 2 2 1 0 1 2 0 2 1 0 1 2 1 2 1 ) /( 1 ) ( ) ( . (D17)
The vector  ~ which minimizes the quantity:

  2 1 ) ( ) (     K k T t t y J   (D18) with a discretization k t , K k ,..., 1  of t, is given by: ) ( ) ( ) ( ) ( ~1 1 1 k K k k k T K k k t y t t t                 . ( D19
)
We then have the following relations: 1 and the multiplicative factor  (=5.015) is equal to the ratio of the distance between two consecutive immobile zones to the radius of gyration of the immobile domain to the mobile zone.  is the ratio of the total immobile porosity to the mobile porosity. for the MRMT models equivalent to the four diffusive porosity structures presented in Figure 2, as determined by the diagonalization method (section 3). Normalized rates larger than 200 or corresponding to a normalized porosity smaller than 10 -3 have been truncated. .

  Figure 1: (a) Skeleton of a dissolution feature in an oolitic limestone, observed by X-ray micro-tomography [Luquot et al., 2014]. The dissolving acidic solution percolates from top to bottom on the general view (bottom left). Its pH increases from top to bottom and from inside out of the main flow path indicated by the curved arrow on the detailed view (top right). The acid dissolves preferentially the calcite cement surrounding the oolites, the size of the pores progressively decreases away from the main flow path, and the organization of the pores becomes more complex. (b) Structured INteracting Continua model (SINC) sketched from the dissolution pattern of (a) with three cross sections transversal to the mobile zone materialized by the arrow. (c) Equivalent MRMT model with the 5 most important rates as determined by the numerical methods set up in section 4. The size of the boxes scales with the porosity affected to the rates labeled by triangles in Figure 7.

Figure 2 :2

 2 Figure 2: Examples of Structured INteracting Continua (SINC) used to illustrate and validate the numerical identification methods of the equivalent MRMT models. From left to right, the diffusive porosity structures are (a) the classical Multiple INteracting Continua (MINC) [Pruess and Narasimhan, 1985], (b) an asymmetric Y with a single junction, (c) an asymmetric loop, and (d) the dissolution structure presented in Figure 1. The size of the immobile cells is proportional to their porosity and the distance along the immobile structure is to scale. The mobile zone is represented by the thick black box with the crossing arrow. Its size has been exaggerated 10 times to be clearly marked. To be comparable, the four structures have the same total porous volume and the same radius of gyration taken with respect to the mobile zone.

Figure 3 :

 3 Figure 3: Diffusive porosity structures represented as cross-sections transversal to the mobile zone direction ((a),(b)), with their associated interaction matrix A ((c),(d)) for the asymmetric Y (top) and MRMT structures (bottom). Dotted frames around subsets of the immobile porosity structures ((a) and (b)) and around matrix lines ((c) and (d)) show how structures are translated in matrix form. Parameters for the asymmetric Y structure are taken from Table 1

Figure 2 .

 2 Figure 2. Because of their widely differing structures, these four SINC examples are thought to be a good basis for testing and illustrating the numerical methods. The first one is the classical MINC taken as reference. The two next ones were chosen for their elementary branching and looping connectivity patterns (Figure 2b and c). Any more involved patterns like the dissolution pattern of Figure 2d will be some kind of combination of these elementary structures. These four examples are comparable in the sense that they have the same mobile

Figure 4 :

 4 Figure 4: Diffusive porosity structures used to check the numerical implementation of the SINC model, equivalent to the "MINC 1D" structure (left column). These structures display the same behavior for homogeneous initial concentrations in the immobile zones. The mobile zone is the bold box with the arrow. The size of the boxes is proportional to the porosity of the compartments. Only the vertical distance of an immobile zone to the mobile zone is to scale.

Figure 6 :

 6 Figure 5: Differences in macrodispersions ) , ( MRMT SINC D D diff(equation (20)) between SINC models and their approximate MRMT models with a limited number of rates 5 to 1 

  corresponding to the immobile zones connected to the mobile zone, due to the removal of the column corresponding to the mobile zone in the extraction of real, having only negative diagonal elements, and also being irreducible diagonally dominant. The diagonally dominance has been shown previously. The irreducibility property is more involved but can be proved by studying the properties of the graph defined by structure connected to the mobile structure by a single link, at least one path exists from any immobile cell to any other immobile cell that does not cross the mobile zone, then the graph defined by M  is strongly connected, so M  is irreducible[START_REF] Horn | Matrix analysis / Roger A. Horn[END_REF]. When several immobile zones are independently connected to the mobile zone, each of these immobile zones is associated to a strongly connected graph and to an irreducible diagonally dominant matrix, itself a sub-matrix of ones of the sub-matrices.

  C2), whose coefficients are given by  1 .. N . The porosities i im  are further deduced from  N+1 .. 2N by inversing the N+1..2N equations of (C7):

  deduced from the identified values of i  . In the case of strict equivalence between MRMT and SINC models, the equivalent MRMT model can be found through (C6)-(C9).In the case where MRMT and SINC models are not strictly equivalent, we seek for the composition of N exponential functions that best matches )

N

  is the largest rate, the initial time sampling should be smaller of Shannon's theorem. Adequate time sampling could then increase with time for determining the smaller rates.

  in the Laplace domain, we separate the time-dependent elements from the quantities depending on i im

  We recall the expression of the discharge of one immobile zone in MRMT model into a mobile zone of constant concentration zero (equation (C1)): t) is the remaining mass per unitary volume of solute and 0 c is the initial homogeneous immobile concentration.

Figure 1 :

 1 Figure captions

Figure 2 :

 2 Figure 2: Examples of Structured INteracting Continua (SINC) used to illustrate and validate the numerical identification methods of the equivalent MRMT models. From left to right, the diffusive porosity structures are (a) the classical Multiple INteracting Continua (MINC) [Pruess and Narasimhan, 1985], (b) an asymmetric Y with a single junction, (c) an asymmetric loop, and (d) the dissolution structure presented in Figure 1. The size of the immobile cells is proportional to their porosity and the distance along the immobile structure is to scale. The mobile zone is represented by the thick black box with the crossing arrow. Its size has been exaggerated 10 times to be clearly marked. To be comparable, the four structures have the same total porous volume and the same radius of gyration taken with respect to the mobile zone.

Figure 3 :

 3 Figure 3: Diffusive porosity structures represented as cross-sections transversal to the mobile zone direction ((a),(b)), with their associated interaction matrix A ((c),(d)) for the asymmetric Y (top) and MRMT structures (bottom). Dotted frames around subsets of the immobile porosity structures ((a) and (b)) and around matrix lines ((c) and (d)) show how structures are translated in matrix form. Parameters for the asymmetric Y structure are taken from Table 1

Figure 4 :

 4 Figure 4: Diffusive porosity structures used to check the numerical implementation of the SINC model, equivalent to the "MINC 1D" structure (left column). These structures display the same behavior for homogeneous initial concentrations in the immobile zones. The mobile zone is the bold box with the arrow. The size of the boxes is proportional to the porosity of the compartments. Only the vertical distance of an immobile zone to the mobile zone is to scale.

Figure

  Figure 5: Differences in macrodispersions ) , ( MRMT SINC D D diff (equation (20)) between SINC models and their approximate MRMT models with a limited number of rates 5 to 1  N , for the four SINC models presented in Figure 2. The determination of the approximate MRMT models is achieved with the numerical identification method in the temporal domain (Section 4.1 and Appendix C).

Figure 6 :

 6 Figure6: Breakthrough curves for the dissolution-like SINC model (Figure2d) and for its equivalent MRMT models, either determined by the diagonalization method (section 3), or by the numerical method in the temporal domain with a limited number of N rates (Section 4.1 and Appendix C). The concentrations are measured at the position

Figure

  Figure 7: Normalized rates   i versus normalized porosities

1

Acknowledgements

The ANR is acknowledged for its funding through its project H2MNO4 under the number ANR-12-MONU-0012-01. The authors are also grateful to Linda Luquot, Alain Rapaport and Jesús Carrera for stimulating discussions, and additionally to Linda Luquot for providing the illustrative dissolution patterns. We also thank the four anonymous reviewers for their insightful and critical review of the manuscript, as well as Cass Miller for his editorial work.

Appendix C: Numerical identification method of MRMT models equivalent to a SINC model

We set up an optimization scheme to get MRMT models equivalent to a SINC model. We first consider the case where the mass per unitary volume m(t) discharged to an immobile zone from a flushing experiment can effectively be modeled by a series of N exponential functions with rates  i and associate porosities

, see equation ( 15)). We derive a set of equivalent expressions in the Laplace domain with simple dependences on  i and i im  . We then deduce optimization strategies both in the Laplace and temporal domains.

We assume first that the SINC model is strictly equivalent to a given MRMT model as in

). It is the case when initial concentrations are homogeneous in the