Solute transport in heterogeneous reservoirs: Upscaling from the Darcy to the reservoir scale

Accepted Manuscript

Solute transport in heterogeneous reservoirs: Upscaling from the Darcy to the reservoir scale Carlos G. Aguilar-Madera , Erik C. Herrera-Hernandez , ´ Gilberto Espinosa-Paredes PII: DOI: Reference:

S0309-1708(18)30634-1 https://doi.org/10.1016/j.advwatres.2018.12.002 ADWR 3254

To appear in:

Advances in Water Resources

Received date: Revised date: Accepted date:

23 July 2018 10 November 2018 12 December 2018

Please cite this article as: Carlos G. Aguilar-Madera , Erik C. Herrera-Hernandez , ´ Gilberto Espinosa-Paredes , Solute transport in heterogeneous reservoirs: Upscaling from the Darcy to the reservoir scale, Advances in Water Resources (2018), doi: https://doi.org/10.1016/j.advwatres.2018.12.002

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights 

We present the averaged mass transport equations for solute in two-facies reservoirs



The two-equation model was developed along its closure scheme and effective coefficients We computed numerically effective coefficients for lenticular and layered reservoirs

AC

CE

PT

ED

M

AN US

CR IP T



1

ACCEPTED MANUSCRIPT

Solute transport in heterogeneous reservoirs: Upscaling from the Darcy to the reservoir scale Carlos G. Aguilar-Madera*a, Erik C. Herrera-Hernándezb, Gilberto Espinosa-Paredesc a

CR IP T

Universidad Autónoma de Nuevo León, Facultad de Ciencias de la Tierra, Ex-Hacienda de Guadalupe, C.P. 67700, Linares, N.L. México. Corresponding autor: [email protected] b

CONACYT-Centro de Ingeniería y Desarrollo Industrial, Av. Playa Pie de la Cuesta 702, Desarrollo San Pablo 76125, Querétaro, Qro., México

División de Ciencias Básicas e Ingeniería, Universidad Autónoma Metropolitana-Iztapalapa, C.P. 09340, Del. Iztapalapa, Ciudad de México, México

AN US

c

Abstract

In this work, we present a methodology for upscaling solute transport in heterogeneous reservoirs. The problem considers upscaling from the Darcy scale of a porous medium to the reservoir scale. We first define the mathematical problem that governs solute mass

M

transport at the so-called Darcy scale, and then we apply the method of volume averaging. This method allows for upscaling the governing equations at the whole reservoir scale and,

ED

eventually, developing physical constraints, definitions of average variables and effective coefficients. One important characteristic of the volume averaging method is that the so-

PT

called closure problems for analytical or numerical computation of effective coefficients are identified. The closure problems consist of boundary-value problems, whose solution

CE

must be carried out in representative geometries (unit cells) of the reservoir. In this paper, we developed the theory for a two-facies reservoir, but the extension for more facies is straightforward. In this context, we present the two-equation model, and the numerical

AC

estimation of effective coefficients are presented and discussed for lenticular and layered arrangements of facies. In this manner, relevant coefficients are identified. The numerical solution of the upscaled model agrees well with direct numerical results when the governing equations are solved at the Darcy scale. Keywords: mass dispersion; reservoir scale; volume averaging method; two facies, twoequation model 2

ACCEPTED MANUSCRIPT

1. Introduction The modeling of solute mass transport in heterogeneous reservoirs has been an active theme for investigation. This theme is of interest in areas such as hydrogeology, enhanced hydrocarbon recovery, environmental sciences, as well as geothermal energy. Commonly, classic advection-dispersion equations have been extensively used to model the transport of

CR IP T

solutes through heterogeneous reservoirs, but the physical meaning and rigorous link of effective coefficients with physics at smaller scales are unclear.

Quintard and Whitaker [1] considered the transport of an adsorbing solute in a two-region model of a chemically and mechanically heterogeneous porous medium when the condition

AN US

of large-scale (for instance, at the field scale) mechanical equilibrium is valid. The authors consider that under these circumstances, a one-equation model can be used to predict the large-scale averaged velocity, but a two-equation model may be required to predict the regional velocities that are needed to accurately describe the solute transport process. If the condition of large-scale mass equilibrium is valid, the solute transport process can be

M

represented in terms of a one-equation model. In this context, on one hand, the twoequation model refers to a set of governing equations composed of one equation to model

ED

flow inside one region, and the other equation models fluid flow in the remaining region. In this formulation, the properties of each region are required. On the other hand, for the oneequation formulation, only one equation is employed that contains global parameters, and it

PT

is not able to model fluid flow in individual regions. The upscaling procedure followed by Quintard and Whitaker consisted of scaling the point equations (i.e., local-scale) to the

CE

Darcy scale (i.e., to the scale of the porous medium); then the authors scaled the Darcy description to a larger-scale, applying the volume average method twice, which represented

AC

a method of systematic upscaling based on upscaling the governing equations in sequential form. On the other hand, Porta et al. [2] presented a volume-averaging-based theoretical analysis of upscaling of reactive transport processes involving fast, bimolecular, homogeneous, irreversible reactions occurring within a porous medium. The authors began from the formulation that drove system dynamics at the pore scale and derived the governing equations at the observation scale associated with laboratory-scale scenarios.

3

ACCEPTED MANUSCRIPT

This involved advection-dominated transport of two reactants and resulting products in the presence of different relative strengths (Dahmköhler and Péclet numbers). In the seminal work by Wood [3], a valuable discussion is presented on the role of scaling laws in upscaling. He discussed the process of coarse-graining in complex subsurface hydrologic systems with the intention to examine the difference between the mathematical

CR IP T

process of averaging and the upscaling process. According to this author, an averaged description does not contain less information than the original problem. For the problem of nonreactive solute transport in heterogeneous media, five different scaling laws were identified as necessary: (i) influence of boundaries, (ii) statistical homogeneity, (iii) separation of length scales, (iv) ergodicity, and (v) smallness of the variance. Another

AN US

discussion by Wood concerns an upscaled mathematical model that has the ability to represent a complex system at the expense of also requiring more parameterization and characterization in order to get the model to work. Regarding the scales, porous systems often exhibit hierarchical configurations [e.g., 4,5], where for instance, the solid matrix at the macroscopic scale is itself microporous. Another example concerns porous media in

M

which the impermeable solid mesostructure is embedded in a finer, saturated, porous medium, which is also embedded in another finer porous structure, and so on. For this

ED

reason, heterogeneous reservoirs can be considered at the following hierarchy of scales: local scale, Darcy scale, macroscopic scale and large scale or reservoir scale [4]. This

systems.

PT

categorization provides a scheme for applications in geo-energetic and hydrogeological

CE

Raoof and Hassanizadeh [6] developed relationships between the pore-scale adsorption coefficient and corresponding upscaled adsorption parameters. They applied both

AC

theoretical averaging and numerical upscaling. In the averaging approach, equilibrium adsorption is assumed at the pore-scale, and solute transport equations are averaged over a representative volume (REV). This leads to explicit expressions for macroscale adsorption rate constants as a function of microscale parameters. In the numerical approach, Raoof and Hassanizadeh first simulated solute transport within a single tube undergoing equilibrium adsorption at the pore wall, from which flux averaged concentration breakthrough curves

4

ACCEPTED MANUSCRIPT

were obtained. These were used to determine the upscaled adsorption rate constants as functions of the pore-scale. Other upscaling methods for solute transport have been proposed. For example, a method based on genetic programming was presented by Hill et al. [7], and a method based on random walking was proposed by Geiger et al. [8]. Hill et al. investigated the development

CR IP T

of a domain-independent modeling tool that searches the space of mathematical equations for one or more equations to describe a set of training data. These methods, based on genetic algorithms, are promising when there are many field data, but in most cases when this information is available, the method is no longer necessary; that is, this method can be applied for post-analysis for interpretation of physical properties. Currently, classical dual‐

AN US

porosity approaches in conjunction with an advection‐dispersion equation and macroscopic dispersivity commonly fail to accurately predict the breakthrough of fractured porous media; for this, a continuous time random walk method can be used as a generalized upscaling approach [8]. According to Geider et al., this procedure allows for predicting breakthrough at other locations accurately and gaining significant insights into the nature of

M

fracture‐matrix interactions in naturally fractured, porous reservoirs with geologically realistic fracture geometries. However, this method requires the correct estimation of the

solute transport.

ED

distribution of retention times, which is crucial for understanding the upscaled character of

PT

The upscaling of mass dispersion has also been studied using the percolation theory in combination with scaling methodologies. Recently, Sahimi [9] employed the continuous-

CE

time random walk method along the power-law growth of the mean square displacement of the solute to analyze mass dispersion in heterogeneous media. Sahimi presented the

AC

intrinsic relation of geometrical exponents of percolation with the coefficients involved in the power-law of effective conductivity and permeability close the percolation threshold. Ghanbarian-Alavijeh et al. [10] extended the analysis of dispersion involving the saturation and porosity effects. The studies combined the percolation technique with fractal scaling laws. They found that the most important feature for facilitating good experimental predictions is the applicability, or not, of random or invasion percolation simulations.

5

ACCEPTED MANUSCRIPT

In enhanced hydrocarbon recovery studies, combustion tube experiments are conducted in the laboratory to determine the burning characteristics of a particular oil in a rock matrix. The data obtained from history matching the tube experiments cannot be directly used for field simulation because grid blocks for the tube simulation are on the order of inches thick, whereas grid blocks for a field simulation are typically 100 to 1,000 times larger [11]. As

CR IP T

such, an upscaling procedure is necessary. In this sense, on one hand, Marjerrison and Fassihi [11] presented a procedure for scaling heavy-oil combustion tube results to a field model. The authors showed that their scaling method for a simulation model was feasible. However, it is not a systematic method because it is based on empiricisms. On the other hand, Islam and Ali [12] analyzed the scaling criteria for in-situ combustion experiments.

AN US

They used partial differential equations in three-dimensional systems to derive a set of scaling criteria. These criteria can be applied, but they do not constitute an upscaling method. Recently, Aguilar-Madera et al. [13] presented an upscaling procedure for in-situ combustion by using the volume averaging method. Considering homogeneous and heterogeneous chemical reactions, multiphase and multicomponent features, the authors

M

developed closure schemes to estimate effective coefficients, but the application of such theory has not yet been done. On the other hand, the recovery from incompletely water-wet

ED

fractured reservoirs can be extremely low. A reason for low recovery is related to wetting issues, whereas the reason for slow recovery may be the nonequilibrium behavior of capillary pressure. To incorporate nonequilibrium into larger-scale problems, Salimi and

PT

Bruining [14] applied homogenization techniques to derive an upscaled model for fractured reservoirs. The association of methods for the complex scaling problem may be a feasible

CE

solution; for example, the association of upscaling by the volume averaging method with the network approach is presented in the work of Ahmadi et al. [15].

AC

The facies analysis approach is widely used and accepted to study heterogeneous reservoirs, as well as in fields such as deep-water processes [e.g., 16-18], glacier formation [e.g., 19-21], diagenesis in oil fields [e.g., 22-23], carbonate sedimentary systems [e.g., 2425], other lithological studies [e.g., 26], pore characterization [e.g., 27], and heterogeneity and compartmentalization within a reservoir for seismic facies analysis [e.g., 28-29]. In the work of Cullis et al. [18] facies associations were applied to diagnostic criteria used to discern each hierarchical order in the field of deep-marine sedimentary geology to assign a 6

ACCEPTED MANUSCRIPT

spatial and temporal order to the sedimentary architecture of preserved deep-marine deposits. Meanwhile, in [20], ground penetration radar (GPR) was used to identify sedimentary facies. The detailed heterogeneity of shale facies and their relationship to pore networks was analyzed in Li et al. [24]. Various approaches are used to develop numerical models for facies distributions within a

CR IP T

reservoir. However, only two approaches are traditionally used [30]: pixel-based and object-based techniques. It is worth mentioning that both approaches can be used simultaneously for facies distribution numerical models. In a recent study, advanced modeling approaches for gas field analysis considered multistage constraints and

AN US

hierarchical facies control as well as multistep modeling [31].

According to previous studies, facies analysis is fundamental because porosity and permeability control the storage and fluid flow in reservoir rocks, which are highly heterogeneous due to their geological evolution, influenced by sedimentary origin, diagenetic processes, and burial history. The latter processes influence the size and shape of

M

pores that produce complex pore networks in sedimentary rocks. In this regarding, lattice Boltzmann simulations of fluid flow in continental carbonate reservoir rocks was presented

ED

by Soete [32], where dominant lithofacies types were identified. All of the facies types had pores that varied from at least the nanometer to centimeter scale, i.e., meso- and macropore networks, which had relevant contributions to fluid flow in the material. Thus, facies

PT

modeling is an essential part of reservoir characterization where the connectivity of the facies model is critical for the dynamic modeling of reservoirs. According to Bayat [33],

CE

carbonate reservoirs are so heterogeneous that variogram-based methods such as sequential indicator simulations are not very useful for facies modeling. The author also proposes the

AC

use of multiple point geostatistics (MPS) for facies modeling in oil fields. In this work, we propose a new methodology for modeling solute transport in heterogeneous reservoirs composed by two-well defined facies from the point of view of own petrophysical properties. The fundamental idea behind this methodology is the upscaling of the model, where the Darcy scale corresponds to porous media and the application of a volume averaging procedure allows upscaling of the governing equations for solute transport at the reservoir scale. The upscaling method generates more unknown 7

ACCEPTED MANUSCRIPT

variables than equations, which are solved with closure problems in representative geometries to obtain the effective coefficients at the reservoir scale. A facies analysis of the reservoir must already exist, i.e., the starting point of the methodology developed in this work considers a known hierarchical facies structure. Although the theory presented in this work corresponds to a two-facies reservoir, the extension to a larger number of facies could

CR IP T

be also carried out. However, the application to three or more facies reservoirs is extremely complex, and it will yield three or more averaged equations strongly coupled through several effective coefficients. The application in practical problems of such mathematical model would be an overwhelming task, from the point of view of the model solutions and the determination of the various effective coefficients involved. Nevertheless, we state that

AN US

a simple one-equation model can be obtained adding the individual averaged equations for each facies, and this model could be more representative for three or more facies reservoirs. The extension of the volume averaged methodology here presented to three or more facies is out of the pursued scopes, and it will be investigated in future works. For the two-facies approximation, the numerical estimation of effective coefficients is presented and discussed

M

for lenticular and layered arrangements of facies. According to results obtained with a numerical solution, the upscaled model agrees well with direct numerical results when

ED

governing equations are solved at the Darcy scale. This article is organized as follows: Section 2 presents the governing equation, at the Darcy

PT

scale, for two well-defined facies with homogeneous properties in the reservoir. Section 3 describes the preliminary results of the volume averaging procedure. In Section 4, the

CE

closed system of equations in terms of the spatial deviations of concentrations is developed. In Section 5, the closed volume averaged equations and the physical meanings of effective coefficients are discussed. Section 6 presents the numerical estimation and analysis of

AC

effective coefficients. In Section 7, the upscaling model validation (reservoir scale) is performed through comparison with the Darcy scale simulation. Major observations regarding the proposed methodology for upscaling modeling (from the porous scale to the reservoir scale) are presented in Section 8.

2. Governing equations at the Darcy scale 8

ACCEPTED MANUSCRIPT

We begin the study of solute transport at the Darcy scale by assuming that there are two well-defined facies with homogeneous properties in the reservoir. As shown in Figure 1, for one heterogeneous reservoir composed of two homogeneous facies, called the  - and  regions, the accommodation of porous rocks inside the reservoir can be one of two categories: (a) dispersed, and (b) layered two-facies reservoir. The origin of these

CR IP T

geological types of reservoirs is driven by, for instance, the sedimentary environment and sediment provenance. For layered reservoirs, only one configuration can be found, but for dispersed two-facies reservoirs, there are four variants, namely, (a) flaser, (b) wavy, (c)

AC

CE

PT

ED

M

AN US

lenticular and (d) ―pin stripe‖ types, as depicted in Figure 1.

Figure 1. Scheme of one heterogeneous reservoir with two well-defined facies (regions). Here, the characteristic lengths of volume averaging and homogeneous regions are depicted. 9

ACCEPTED MANUSCRIPT

For the two-region reservoirs depicted in Figure 1, the system of equations governing the mass transport of one solute inside a mobile fluid phase is given by: In the  -region

CR IP T

C     u C      D C  t In the  -region C







   u C    D C



AN US

t

At t  0

C  C 0 ,

At the    boundary

C  C 0

M ED

(2)

(3)

C  C

At the    boundary

(1)

(4)





(5)

PT

n   u C  D C   n  u C  D C

In these equations Ck and Dk are the solute concentration and mass dispersion tensor in

CE

the k -region, respectively. In the transport problem represented by Eqs. (1)-(5), it is assumed that the seepage velocities u and u can be calculated alternatively according to

AC

the type and number of fluids flowing in the reservoir, as well as with previous knowledge of petrophysical properties of the reservoir. For instance, a Darcy’s Law assumption can be used for liquids (oil and water), but the Darcy-Forchheimer model is commonly employed for gases. It is worth mentioning that dispersion coefficients are intrinsic functions of microstructure and pore space characteristics [34] and fluid velocity [35]. In Eq. (5), n is the unit normal vector pointing from the  -region to the  -region (see Figure 1).

10

ACCEPTED MANUSCRIPT

The boundary-value problem for C is completed by the statement of conditions at external reservoir boundaries. These conditions can be set up as known flux ( J ,re ) or concentration ( C ,re ), which can be position-dependent, i.e.,: At the external boundary of the reservoir (6)

C  C ,re  z 

(7)

CR IP T

n ,re   u C  D C   J ,re  z 

or

AN US

Here, n ,re is the normal vector pointing outside the reservoir and z is the position vector locating the external boundary. Eqs. (6) and (7) are provided for the sake of completeness, but as it is developed further, this information will be not required because boundary effects are discarded in our analysis and representative zones of reservoir bulk are used to predict

M

numerically effective coefficients.

ED

3. Preliminaries: volume averaging procedure In the volume averaging method it is necessary to define superficial and intrinsic volume

PT

averaging operators as follows:

CE

Superficial average

AC

Intrinsic average

k 

1  k dV , V V

k   ,

(8)

k

k

k



1  k dV , Vk V

k   ,

(9)

k

where  k is any variable defined in the k -region, Vk represents the geometrical space only encompassing the k -region where the averaging is carried out, and V and Vk are the volumes of geometrical spaces V and Vk , respectively. The averaging operators (8) and (9) are related by the volume fraction of k -region as follows: 11

ACCEPTED MANUSCRIPT

 k  k  k

k 

k

, k  ,

(10)

Vk , k   , V

(11)

The size and shape of the averaging volume must be appropriate in order to spatially smooth the variable that is to be upscaled. In Figure 1, it is apparent that for dispersed

CR IP T

reservoirs, one spherical shape encompassing various ―  -type rocks‖ might be the best option; meanwhile, for layered reservoirs, one thin rectangular prism including the entire net pay of the reservoir could be the most appropriate.

We first conduct the theoretical development to upscale the solute transport inside the  -

AN US

region. Thus, applying the operator (8) to Eq. (1), we have

C     u C      D C  t

(12)

To introduce the integral operator within the time and space derivatives, we appeal to the

M

spatial averaging theorem [36],

ED

      

1 V



n  dA

(13)

A

PT

and the Reynolds general transport theorem 1   k w   n dA  V   k we   nedA

A

(14)

Ae

CE

    1   t t V

In the last equation, A is the interfacial area contained in V , w and we are

AC

velocities of   boundary and of the external surface ( Ae ) of representative elementary volume (REV) in contact with the  -region. Here REV is equivalent to the averaging volume. ne is the unit normal vector pointing from the  -region outward from the REV, as shown in Figure 2.

12

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 2. Scheme of REV showing normal and position vectors. Considering nondeformable and static reservoir rocks and REV, we have (15)

M

w , we  0

It is important to note that the assumption w  0 implies that mass fluxes at the  

ED

boundary are equal, i.e.,  u  n   u  n (here,  and  are the  - and  facies densities, respectively). This means that there is one impenetrability condition

PT

between the  - and  - regions in the sense that one facies does not encroach on the other.

 C C  t t

(16)

AC

CE

Thus, Eq. (14) simplifies to

Turning our attention to the advective term of Eq. (12) and after using the spatial averaging theorem, we have

   u C     u C 

1 V



n   u C  dA

(17)

A

In the same way, the dispersive term becomes 13

ACCEPTED MANUSCRIPT

   D C 

  1     D    C    V   

1 V



  n C dA  A  

(18)

n   D C  dA

A

CR IP T

Note that Eqs. (17) and (18) contain two types of concentrations: the Darcy scale, C , and the average concentration, C ; then, the spatial decomposition [37]

  



 

(19)

and

space

as

it

is

  t, x    t, x   

directly 

t, x .

AN US

is introduced in the first integral of Eq. (18). This deviation variable   is function of time related

to

the

Darcy

scale

variable,

Thus, after some rearrangements, Eq. (12) can be

written as

   1     u C     D     C   n C dA      t V A     1 1     n  u C  D  C dA   n  u C  D C dA V A V A 





M

 C



ED



i.e.,



(20)





CE

PT

Here, the intrinsic average concentration was employed, as was the equality

  

1 V



n dA

(21)

A

AC

which is a direct result of the spatial averaging theorem (13) when    Constant . Substituting the spatial deviations for concentration and velocity into the advective term of Eq. (20), we have u C   u



C



 u C

(22)

14

ACCEPTED MANUSCRIPT

Arriving to this equation, we discarded the terms containing average of spatial deviations, i.e.

  0

(23)

Such an approximation is valid if

L,

r02

LC1L

(24)

CR IP T

r0

where r0 represents the size of volume averaging, L is the characteristic length of the reservoir and LC1 is the characteristic length over which significant concentration gradient changes occur. The derivation of these inequalities is given in Appendix A.





     u

t 

1 V



 n   u

C C

A





   1     u C    D    C   n C dA      V A     (25)  1   D  C dA   n  u C  D C dA V A

M



 C

AN US

Thus, the volume averaged equation for the  -region is











ED

This equation is impractical because it contains the average concentration and its spatial deviation, and additionally, the surface integrals confer nonlocal characteristics. Thus, the

PT

following stage within the volume averaging method refers to developing the governing equation for spatial deviations. We note that a similar procedure is performed to obtain the

CE

volume averaged equation for  -facies, but this development is omitted here for brevity.

AC

4. The closure scheme To develop a closed system of equations, we must find the governing equation for C . A simple and direct procedure is subtracting Eq. (25) from (1), which leads to:

15

ACCEPTED MANUSCRIPT

C     u C  u t 



C



 u C





    D C    

   1 1     D  n C dA    n  u C    V   V   A A      1  n  u C  D C dA V A



 D  C



 dA

(26)





CR IP T







Here, a homogeneous behavior for   was considered, i.e.,   0 , which means that the volume fraction of  -facies within the averaging volume remains constant as the centroid of the averaging volume moves into the reservoir. It is a somewhat strong assumption

AN US

which is straightforwardly satisfied for periodic arrangements of facies. Nevertheless, in practice, geological reservoirs can be highly heterogeneous within a given depositional environment. Thus, the validity of    0 can be applied only to reservoirs that are statistically homogeneous with no strong changes in the accommodation of each lithofacies, such as in the bulk of geological formations and sufficiently far from the tops, faults, base,

M

and pinch outs of reservoirs. Within the geostatistics framework, there is a type of reservoir that is considered as stationary when stochastic processes are performed to model the

ED

spatial distribution of properties (Kriging procedures). In this context, modeling a stationary reservoir implicates that the mean values and variance of petrophysical properties are constants inside the reservoir. Such implications help qualitatively to satisfy

PT

the requirement of a homogeneous facies distribution. However, a more rigorous analysis

CE

of physical implications behind the expression    0 is beyond the scope of this work. We recognize the importance of exploring this topic in a future study, for which we will

AC

study coupling between geostatistical data and the upscaling of geological reservoirs. Using the spatial decomposition for concentration and velocity, the convective term of Eq. (26) can be rewritten as

u C  u



C



 u C



 u C  u C



 u C



(27)

which leads to

16

ACCEPTED MANUSCRIPT

C     u C  u C t  

1 V







 u C

n  u C





      D C     D   1       V   



 D  C





 dA  V1



 A 

A



  n C dA  A   (28)



n  u C  D C dA

CR IP T

Now, based on the analysis of orders of magnitude and the length-scale constraints presented in the Appendix A, the following inequalities appear and the associated expressions in terms of magnitude estimates are also written: C t

D C O  2   L  

C  O  *   t 

M

And,

(29)

  n C dA  A 

  



  D C

(30)



(31)

ED

  1    D     V  

AC

CE

PT

if:

And,



AN US

if:



  D C

 D C A ,REV  O    V   L

  u C



D C O  2   L  



  u C



  

(32)

(33)

if:

u C  O     L 

u C  O     L 

(34)

And, 17

ACCEPTED MANUSCRIPT







n  D  C

 dA





A



n  D C dA

(35)

A

if

D C  O     L 

(36)

CR IP T

D C  O     L 

The definition of each estimate is provided in Appendix A. Since the left-hand terms of the above equations are negligible compared to the right-hand terms, Eq. (28) can be simplified. Then, convective terms including the spatial deviations of velocity can be written as: 

u





AN US



  u C



 C

 C

  u

(37)

A similar procedure to obtain the spatial deviations C can be applied to obtain the governing equation for C . Thus, the complete boundary-value problems for spatial

M

deviations of solute concentrations in the  - and  -regions are:



ED

In the  -region



  u C  u  C



 C

PT

Source term





  u    D C



Source term





n  u C  D C dA 

C V

A

CE

1  V







n  u dA

(38)

A Source term

AC

In the  -region





  u C  u  C Source term

1  V

 A





 C





  u    D C



Source term



n  u C  D C dA 

C V





n  u dA

(39)

A Source term

18

ACCEPTED MANUSCRIPT

At the   boundary

C    C 



 C



C   

(40)

Source term







CR IP T



n  u C  D C  n  u C  D C

     n  u C  u C    Source term Source term 

(41)

AN US

     n  D  C  D  C   Source term   Source term

The boundary conditions (40) and (41) were derived using spatial decompositions in Eqs. (4) and (5); additionally, the continuity of velocity was considered at the   boundary,

M

i.e.

ED

At the   boundary

u  u

(42)

PT

The boundary-value problem (38)-(41) is a set of linear equations in which several source terms are identified. In the domain where spatial deviations exist, these source terms

CE

behave as constant. These features allow us to formulate a linear combination of particular solutions for C and C , where each solution includes the effect of just one source. The

AC

complete solution is the linear sum of all the particular solutions. Additionally, the particular solutions can be solved assuming a unitary source. In this way, the complete solution (named the formal solution) is the linear combination of individual solutions multiplied by its associated source, i.e.: C  b  C



 b  C



 d C



 d C



 s  C 



 C



  A (43)   

19

ACCEPTED MANUSCRIPT

C  b  C



 b  C



 d C



 d C



 s  C 



 C



  A (44)   

In these equations, the functions b kp , d kp and sk ( k  , ) represent the particular solutions and are called closure variables. They can be computed by solving their associated boundary-value problems resulting from the substitution of Eqs. (43) and (44)

CR IP T

into Eqs. (38)-(41). The resulting boundary-value problem can be divided into five separate problems each one involving two coupled closure variables. These five problems are called closure problems, and they allow the computation of particular solutions with a unitary source. The five closure problems are summarized in Appendix B.

AN US

The arbitrary functions A and A in Eqs. (43)-(44) yield a boundary-value problem where the solutions are:

A  A  constant

(45)

Such a constant [and the entire formal solutions (43) and (44)] are then substituted into the

M

average equation (25), but all of the integrals involving A (and A in the corresponding average equation for  -facies) vanish based on the prevalent conditions that are assumed

ED

up to this part of manuscript, i.e.: homogeneous facies properties, the average of spatial deviations equals zero [see Eq. (23)], and the assumption of one incompressible fluid

PT

saturating the pores. For this reason, the functions A and A do not generate any closure boundary-value problem nor effective coefficients. The algebra for this analysis is omitted

AC

CE

for brevity.

5. The closed volume averaged equations 5.1 The two-equation model By substitution of formal solutions (43) and (44) into the averaged equations (25) and the corresponding one for C



, we obtain the volume averaged equation for each facies. The

algebra to obtain the corresponding equation for the  -facies is provided in Appendix C, 20

ACCEPTED MANUSCRIPT

and here we only summarize the main results. Thus, the closed volume averaged equation for the  -facies is: 



   u C

t

accumulation



convection



         D  C  D  C   ordinary dispersion crossed dispersion     w  C





 w  C

   C

non-conservative ordinary and crossed pseudo-convection

    v C 







  C

CR IP T



 C

 v C



pesudo-absorption terms

  

(46)

conservative ordinary and crossed pseudo-convection

AN US

           k L  C  C    av k L  C  C     mass interchange  conservative mass source 



  

whereas the corresponding closed averaged equation for  -region is given by:

t

accumulation

    u 

C



       D  C   D  C            crossed dispersion odinary dispersion   

M



ED



 C

convection

AC

CE

PT

 w  C



 w  C



non-conservative ordinary and crossed pseudo-convection

    v C 



 v C



   C



  C



pseudo-absorption terms

(47)

  

conservative ordinary and crossed pseudo-convection

          k L  C  C    av k L  C  C     mass interchange  conservative mass source 



  

In these equations, several mass transport mechanisms are identified (convection, absorption, and interregion mass interchange) and quantified through several effective transport coefficients. These are the dispersive coefficients D , D , D , and D ;

21

ACCEPTED MANUSCRIPT

the pseudo-velocity vectors v , v , v , v , w , w , w , and w ; the pseudo-absorption coefficients   and  ; and the mass interchange coefficients k L , k L and k L . The effective coefficients are defined as:

  n b dA   u b A 

  n b dA   u b A 

(49)

D

 1  D   V 

  n b dA   u b A 

(50)

AN US

D

 1  D   V 

(51)

 1  D   V 

  n d dA   u d A 

(52)

v

 1  D   V 

  n d dA   u d A 

(53)

v

 1  D   V 

  n d dA   u d A 

(54)

v

 1  D   V 

  n d dA   u d A 

(55)

PT

ED

v

CE

 n b dA      u b A 

M

 1 D  D    I   V 

AC

(48)

CR IP T

D

 1  D    I   V 

22

ACCEPTED MANUSCRIPT

av k L  

1 V

1  V





n  u s  D s dA

A

 1  D   V 

  n s dA   u s A 

(58)

n   D   u b  D b   dA

 A







M

A



1 V

PT

  

1  V

  

1 V

1  V



n  u b  D b dA





(60)

n  u b  D b  D dA

ED

1  V





(59)

n  u b  D b dA

A

1 V

AN US

1 V

CR IP T

(57)

w  

CE



(56)

 n s dA      u s A 

1  V

AC

n   u s  D s  dA

 1 k L  D   V 

k L

w 

 A

A

 n  u 1  d   D d  dA

A

 n   u d  D d  dA

(61)

A

 n   u d  D d  dA A

 n  u 1  d   D d  dA

(62)

A

Here, I is the identity tensor.

23

ACCEPTED MANUSCRIPT

5.2 About physics of effective coefficients Following a more in-depth inspection of effective coefficients defined in (48)-(62), some statements regarding their physics can be concluded: 

Each effective coefficient depends on certain closure variable, and each closure

CR IP T

variable captures the effect of a specific unitary source over solute transport at the Darcy scale. In this regard, we quote, for instance, the coefficient D , which upscales the mass transport by dispersive mechanisms inside the  -facies when there is a driving source  C



, i.e., an inter-well gradient concentration, which

is commonly found in tracer transport tests. In this manner, each effective

AN US

coefficient is related to one source as follows: 

o Effects of the source  C



are quantified in D , D and w

M

o Effects of the source  C

are quantified in D , D and w

ED

o Effects of the source C o Effects of the source C





are quantified in v , v and  

are quantified in v , v and 

PT

o Effects of the source  C 



 C



 are quantified in a k , k  v L L and 

The fact that each scale source can be individually quantified is due to the clear separation of scales (between the Darcy and reservoir scales) that is assumed during

AC



CE

k L

theoretical development [length-scale disparities, Eq. (24)], as well as the linearization imposed to some quantities during the estimation of orders of magnitude (recall that truncated Taylor-series expansions were used). The final result of estimations of orders of magnitude is summarized through Eqs. (29)-(35).

24

ACCEPTED MANUSCRIPT



The effect of the geometrical arrangement of facies is quantified through the area integrals that are included in the definitions of effective coefficients.



Some coefficients include hydrodynamic Darcy scale dispersive effects, which lead to additional transport due to flow within facies. This phenomenon is quantified

CR IP T

through the volume integrals u f x (where f stands for any facies and x stands for any closure variable). 

Most coefficients ( D , D , v , v , v , v , av k L , w , w ,   ,  ) vanish as long as one facies disappears; meanwhile, others simplify ( D , D ). As result of the upscaling procedure, two first-order absorption-like terms appeared naturally,   C



and  C

AN US





. To avoid confusion, they are named pseudo-

absorption terms, as they do not truly represent absorption because this phenomenon was not considered in the Darcy scale governing equations. The

M

physical meaning of these terms can be deduced from the definitions of effective coefficients   and  in Eqs. (61) and (62). They represent inter-facies mass flux

ED

and its driving forces are the average concentrations 



and

C



. For

 0 ,   vanishes. Note in Eqs. (46) and (47) that

PT

instance, for cases when C

C

the absorption-like terms have opposite signs, indicating that the mass of solute

CE

disappearing from one facies appears in the other facies.

AC

6. Numerical estimation of effective coefficients In this section, we performed the numerical estimation of all effective coefficients defined in Eqs. (48)-(62). To this end, the closure problems presented in Appendix B are solved in representative geometries at the Darcy scale. After the closure variables solution is available, the calculation of effective coefficients is performed by computing the

25

ACCEPTED MANUSCRIPT

corresponding integral quantity of definitions (48)-(62). Each closure problem allows the

Closure problem 1 solves for D , D and w



Closure problem 2 solves for D , D and w



Closure problem 3 solves for v , v and  



Closure problem 4 solves for v , v and 



Closure problem 5 solves for av k L , k L and k L

AN US



CR IP T

computation of the following effective coefficients:

Recall that each closure problem requires parameters at the Darcy scale. If such information is known a priori, then this information can be used in the numerical computation. In this work, we do not intend to predict coefficients using real and measured quantities from geological fields; rather, we use the approximate values for representative lithologies

M

summarized in Table 1. Meanwhile, the synthetic geometries depicted in Figure 3 will be used as cells for the solution of closure problems. We stress that our theoretical

ED

development is not limited to these simplifications, and eventually, real geometries and data from the Darcy scale can be utilized to estimate the upscaled coefficients. Within the volume averaging method, the use of representative cells assumes that such geometries can

CE

PT

be periodically repeated in all directions, and eventually, the entire reservoir is reproduced.

AC

Table 1. Relevant information for the solution of closure problems and estimation of effective coefficients.

Parameter

Values

Permeability of  -facies

10 mD

Permeability of  -facies

According to the ratio:

K K

26

Porosity of  -facies

8%

Porosity of  -facies

1%

Pressure drop in the cell

10 psi

Mass dispersion coefficient

5m

CR IP T

ACCEPTED MANUSCRIPT

Molecular diffusion coefficient 1×10-5 cm2/s Flowing phase

Water

AN US

The cell pressured drop remained fixed to solve the closure problems. In this manner, the momentum and mass balance for a unique flowing phase (water) is solved, allowing the computation of the velocity field, which is required as input data in closure problems, and allowing for recalculating the mass dispersion coefficients at the Darcy scale for each facies

M

( D and D ). The Darcy scale coefficients are computed according to the expression: f  ,

(63)

ED

D f  DmI  Α f u f ,

where Dm is the molecular diffusion coefficient and Α f is the mass dispersion tensor,

PT

which is assumed to be isotropic and the same for each type of rock, i.e., Α f   I . According to Eq. (63), the molecular diffusion Dm is the same for both facies because the

CE

flowing fluid (water) and solute are the same in each facies. In this work, the closure problems were numerically solved using Comsol Multiphysics,

AC

which is based on the finite element method. The number of domain elements is shown in Figure 3 for each cell. We noted that the number of domain elements are sufficient to obtain mesh-independent results while not compromising computation efforts. Other numerical configurations were maintained by default, which is optimized for handling nonlinear problems and small relative error tolerance. We limited the number of cases under study, and only two representative geometries of two-facies systems are used: i) a layered system and ii) a dispersed system (lenticular arrangement). Note that for layered 27

ACCEPTED MANUSCRIPT

systems, the height cell is fixed to 20 m and the corresponding computations are

ED

M

AN US

CR IP T

representative of thin or ultra-thin bed reservoirs.

PT

Figure 3. Synthetic cells used for estimation of effective coefficients. Here, the grid mesh,

CE

lengths, and volume fraction of  -facies are presented.

AC

6.1 Closure variables field In Figures 4-8, closure variables computed in representative cells are plotted for four relevant situations: for layered and lenticular two-facies arrangements, and for when one facies is more hydraulically conductive than the other. This is quantified through the permeability ratio K K . The longitudinal components of closure variables involved in

  x , are depicted in Figure 4. These variables can

closure problem 1, i.e.,  b  x and b

28

ACCEPTED MANUSCRIPT

present symmetric or anti-symmetric behavior around defined orthogonal planes. For the lenticular case and considering a poorly conductive  -facies, a quasi-antisymmetry of the closure variables appears. The numerical solution of closure variables related to closure

  x and  b  x ] has similar features to that for closure problem 1, but with

problem 2 [ b

opposite values when lenticular cells are used, as depicted in Figure 5. However, notice that

CR IP T

color legends for each figure have different value ranges. In general, the numerical solution of closure problems 1 and 2 suggests that relevant information takes place around the  - 

CE

PT

ED

M

AN US

interface.

  x [m] as a function of permeabilities and

AC

Figure 4. Closure variables  b  x [m] and b

facies distribution.

29

AN US

CR IP T

ACCEPTED MANUSCRIPT

  x [m] and  b  x [m] as a function of permeabilities and

Figure 5. Closure variables b

ED

M

facies distribution.

The closure variables related to closure problem 3 ( d and d ) are deployed in Figure

PT

6. It is noted for layered cells, that closure variables practically have negligible values, whereas for lenticular arrangements, significant values and symmetry and antisymmetry

CE

appear depending on the permeability ratio. Considering the closure variables of problem 4 ( d and d ), the general trends of the numerical solution are similar to those of closure problem 3 for a lenticular arrangement, as is apparent in Figure 7. The closure variables for

AC

problems 3 and 4 take major relevance when one facies is dispersed into the other. This suggests that the effective coefficients related to the closure problems could take importance for lenticular reservoirs. By recalling the effective coefficient definitions, we note that closure problems 3 and 4 allow for the computation of pseudo-convective coefficients and pseudo-absorptive terms of up-scaled mass balance equations [given by Eqs. (46) and (47)], where some have negligible effects over the global mass transport process. This topic is revisited later in the manuscript. Finally, the numerical solution of 30

ACCEPTED MANUSCRIPT

closure variables related to closure problem 5 ( s and s ) is presented in Figure 8. In this case, for layered and lenticular cells, the closure variables have non-negligible values regardless of the permeability ratio. As in the other closure variables, some symmetry

ED

M

AN US

CR IP T

planes can be defined.

PT

Figure 6. Closure variables d (dimensionless) and d (dimensionless) as a function of

AC

CE

permeabilities and facies distribution.

31

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 7. Closure variables d (dimensionless) and d (dimensionless) as a function of

AC

CE

PT

ED

M

permeabilities and facies distribution.

32

ACCEPTED MANUSCRIPT

Figure 8. Closure variables s (dimensionless) and s (dimensionless) as a function of permeabilities and facies distribution.

CR IP T

6.2 Numerical estimation of effective coefficients After solving all of the closure problems, the various effective coefficients can be computed. We only compute the longitudinal component of coefficients ( ) xx , if applicable, for two arrangements of facies. The effects of the two-facies distributions that are depicted in Figure 1 will be explored in detail in a future study. In the following

AN US

analysis, the dominant parameters are identified, allowing for the reduction of model parameters and optimizing the analysis of effective coefficients.

In Tables 3-6, the numerical estimations of coefficients are tabulated and compared. Several of the coefficients are normalized with Darcy scale mass dispersion coefficients. Note in the first column that the permeability ratio was changed by one order of magnitude.

M

On the basis of their estimations, we conclude that dominant coefficients are those in Table

ED

2.

As complementary information, Table 6 summarizes the associated Péclet number for each case simulated and for lenticular and layered arrangements of facies. The Péclet number

CE

PT

was calculated as:

Péclet number 

Length cell V u dV Vcell cell Dm

,

u  u or u , depending on position (64)

AC

where Vcell is the volume of the cell and u is the norm of the velocity vector. In general, the Péclet number varies exponentially with the permeability ratio, and with the arrangement of facies in the cell.

33

ACCEPTED MANUSCRIPT

Table 2. Dominant effective coefficients according to permeability and the distribution of facies. Lenticular system

Layered system

K  K

 D  xx ,  D  xx

 D  xx

 D  xx ,  D  xx

 D  xx

CR IP T

Facies permeability

 v  x ,  v  x

  ,  K  K

AN US

 v  x ,  v  x  D  xx ,  D  xx

 D  xx

 D  xx

 D  xx

M

 v  x

ED

  , 

PT

av k L

CE

Table 3. Numerical prediction of effective dispersion tensors as a function of permeability

 D  xx  D xx

AC

K K

ratio and facies distribution.

D  D   xx

 xx

D  D   xx

 xx

Layered

D  D   xx

Lenticular

 xx

Lenticular

Layered

Lenticular

Layered

Lenticular

Layered

1×10-4

0.07011

0.50

0.87486

0.50

0.02697

7.98×10-12

0.02808

1.21×10-7

1×10-3

0.06542

0.50

0.87749

0.50

0.02437

9.42×10-12

0.03316

1.24×10-7

1×10-2

0.05780

0.50

0.88678

0.5

0.01753

2.40×10-11

0.04550

1.33×10-7

34

ACCEPTED MANUSCRIPT

0.05801

0.50

0.88807

0.50

0.01664

1.54×10-10

0.04668

1.27×10-7

1

0.06985

0.50

0.87658

0.50

0.02524

7.20×10-10

0.02832

7.20×10-8

1×101

0.09178

0.50

0.86466

0.50

0.03947

1.29×10-9

0.00658

1.49×10-8

1×102

0.10313

0.50

0.65037

0.50

0.18893

1.42×10-9

0.00738

1.68×10-9

1×103

0.09812

0.50

0.85881

0.50

0.03645

1.44×10-9

7.32×10-5

1.71×10-

0.09818

0.50

0.85877

0.50

1.44×10-9

0.03943

AN US

1×104

CR IP T

1×10-1

7.33×10-6

10

1.71×1011

Table 4. Numerical prediction of conservative convective-like coefficients as a function of permeability ratio and facies distribution.

v   D 

 v x  D xx

 x

(1/m)

Lenticular

Layered

Lenticular

3.19×10-4

1.58×10-

3.21×10-4

5.41×10-

3.21×10-4

5.43×10-

3.20×10-4

0.00242

1.56×10-

0.00249

PT

0.00653

1.44×10-

0.00746

10

AC

CE

0.00743 0.00619 0.00186

5.02×10-

0.00864 0.00803

5.95×10-

0.00153

0.00019

8.24×1011

0.00643

9.38×102.35×107.22×10-

2.30×10-4

7.31×10-

0.00721

5.92×1012

9.08×107.35×10-

0.00803

2.35×10-

0.00752

8.22×10-

0.00880

8.66×10-

0.00619

6.90×1012

9.12×102.35×1010

0.00206

8.81×1011

0.01136

12

1.01×10-4

2.11×10-

11

11

3.47×10-4

2.14×10-

10

10

0.00144

2.12×10-

10

11

12

1.30×10-4

0.00246

12

11

11

1×103

1.56×10-

1.34×10-

Layered 10

12

10

10

0.00794

0.00246

12

11

10

1.32×10-

1.57×10-

13

11

11

2.35×10-

13

 xx

(1/m)

Layered

1×10-2

1×102

(1/m)

Lenticular

10

1×101

 x

Layered

1×10-3

1

v   D 

 xx

Lenticular

10

1×10-1

 x

(1/m)

ED

1×10-4

v   D 

 xx

M

K K

1.08×1011

6.67×10-4

1.48×1012

35

ACCEPTED MANUSCRIPT

1×104

4.05×10-

0.00265

1.04×10-4

2.05×10-

10

8.97×10-5

2.49×10-

11

6.55×10-4

5.30×10-

11

13

Table 5. Numerical prediction of non-conservative convective-like and absorption coefficients as a function of permeability ratio and facies distribution.

(1/m) Lenticular 1×10-4

w  D   x

   D  xx (1/m2)

 xx

(1/m)

Layered

4.98×10-4

0.05

Lenticular

Layered

5.79×10-4

4.97×10-

Lenticular

0.00320

0.05

0.00403

4.97×10-

3.94×10-

5.79×10-5

4.54×10-

5.70×10-4

0.05

0.00779

5.02×10-

0.00125

0.05

0.00711

4.68×10-

0.00137

0.04997

0.00316

2.49×10-

1×102

1.37591

1×103

0.04995 0.04995

PT

0.04574

1×104

M

0.04995

0.02995

ED

0.04403

0.04995

0.04333 0.04342

2.20×10-6

0.00306

2.44×10-6

0.04401

7.73×10-7

0.00178

2.62×10-8

6.70×10-4

3.64×10-

8

7

4.33×10-

1.97594

6

9.84×10-

1.11428

9

4.73×10-

0.00315

8

6

8

4.72×10-

1.12×10-6

3.15×10-

0.02290

7

4.69×10-

4.29×10-4

10

1.22×10-4

8.96×10-

6

6.42×10-

11.0787

10

12

9.74×10-5

6.13×10-

5

13

CE

0.04585

0.57486

4.34×10-

8.28×10-7

7.78×10-

0.00465

6

1×101

6.92×10-6

2.91×10-

8.63×10-4

6

1

Layered

1.02×10-

0.00111

6

1×10-1

Lenticular

9

AN US

0.00287

 xx

9

6

1×10-2

D 

(1/m2)

Layered

6

1×10-3



CR IP T

 w  x  D xx

K K

AC

Table 6. Numerical prediction of inter-facies mass transfer coefficients and Péclet number

K K

as a function of permeability ratio and facies distribution.

av kL

D 

 xx

(1/m2) Lenticular 1×10-4

6.02×10-5

Layered 0.03128

kL  D  xx (1/m) Lenticular 3.19×10-4

Layered 1.46×10-

kL

D 

 xx

(1/m) Lenticular 3.20×10-4

Layered 4.83×10-

Péclet number Lenticular 0.09492

Layered 340.264

36

ACCEPTED MANUSCRIPT

7

1×10-3

7.69×10-4

0.03113

0.00242

1.38×10-

6

0.00246

4.84×10-

7

1×10-2

0.00265

0.03023

0.00653

1.29×10-

0.00752

4.82×10-

7

1×10-1

0.00221

0.02675

0.00743

4.63×10-

0.01417

0.00568

2.46×10-

0.00880

0.00125

0.00254

0.00130

4.55×10-

0.00568

4.52×10-

0.00873

2.75×10-

0.00227

4.97×10-

4

0.00162

2.78×10-

0.00133

1.87×10-4

1×104

0.00162

2.78×10-

5.02×10-

0.00332

0.00265

5.02×106

680.460

680.460

4.36×10-

6043.02

3742.53

1.84×10-

59452.0

34363.2

5.94×105

3.40×105

5.93×106

3.40×106

8

6.66×10-4

6

5

374.253

7

6

5

2.52×10-

2.76×10-

AN US

1×103

86.7038

6

6

1×102

343.632

6

6

1×101

9.38922

CR IP T

0.00289

340.570

6

7

1

0.94820

6

6.55×10-4

8

3.22×108

M

Finally, with the goal to observe the general trends of effective mass dispersion and transfer coefficients, they were plotted in Figure 9 as function of the permeability ratio. From the

coefficient

 D  xx .

ED

resulting curves, it is clear that all dispersion coefficients are comparable except the crossed Interestingly, a weak dependency on the permeability ratio is

PT

observed. Regarding the mass transfer coefficient, it is enhanced as long as the permeability

AC

CE

of continuous facies increases, and it is more obvious for lenticular geometries.

37

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 9. Comparison of effective coefficients for layered and lenticular cells, and as a

ED

function of permeability ratio.

PT

7. Validation of model and comparison of dispersion coefficients

CE

7.1 Comparison of model results with Darcy scale simulation The derived two-equation model, Eqs. (46)-(47), is not commonly used to interpret tracer

AC

injection tests in geological reservoirs. To the best of our knowledge, there are not appropriate data to compare and validate numerical estimations of effective coefficients. Instead, the one-equation model is widely used, which only contains one mass dispersion tensor and a unique average concentration. From the two-equation model, it is possible to derive the one-equation model if validity of local mass equilibrium is assumed. However, the implications of local mass equilibrium are unclear at this point; this will be investigated in future work. 38

ACCEPTED MANUSCRIPT

One manner to verify the congruency of the upscaling procedure and whether the effective coefficients were correctly computed is by comparing the solution provided by the upscaled model with those from the Darcy scale governing equations for each facies, Eqs. (1)-(5). We must recall that we initiated the analysis by setting governing equations at the Darcy scale, and then the volume averaging method was applied to upscale the equations. Thus,

CR IP T

we expect similar solutions if the Darcy scale or upscaled models are employed. Nevertheless, it should be clarified that identical solutions are not expected, as some approximations were assumed; i.e., we applied length-scale constraints and estimations of orders of magnitude, as described in Appendix A. Our intention is to demonstrate that the upscaled model has high congruency with the physics that occur at a smaller scale. The

AN US

usage of Darcy scale models involves an overwhelming computational task, which in most of the cases is unnecessary and perhaps unsolvable. As demonstrated later, the upscaled two-equation model has important computational advantages as long as effective coefficients are available.

With these premises, we simulate 30 years of a mass transport case in a 2D, 200×100 m

M

domain, as depicted in Figure 10. We select a lenticular synthetic geometry composed of two-facies with the following petrophysical properties: K  10 mD , K  0.1 mD ,

ED

  0.08 , and   0.01 . Figure 10 corresponds to one transverse cutting plane where the vertical axis represents reservoir depth. At the inlet boundary, solute is injected for a short

AC

CE

as follows:

PT

period at a constant pressure; the outlet boundary is also maintained at a constant pressure

2000 mg/L, t  2 h Csolute   0 mg/L, t  2 h Pinjection  5000 psi Pproduction  4900 psi

39

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 10. Domains for the solution of Darcy scale governing equations, Eqs. (1)-(5), and the upscaled two-equation model, Eqs. (46)-(47). Here, the location of observation points

M

and the grid mesh are also included.

Both Darcy scale and upscaled mathematical models were solved with Comsol

ED

Multiphysics using the grid mesh shown in Figure 10. For a Core i7 2.7 GHz laptop computer, the computational time for the Darcy scale model was 93 s, whereas the upscaled

PT

model required 6 s (93.5% less computational effort). This supports the computational advantages of using average models.

CE

The model solutions were compared by observing the breakthrough curves at observation points identified in Figure 10, as well as the evolution of concentration profiles in

AC

continuous facies (see Figure 11). At first glance, close similitude between evolving concentration profiles is met, but some differences appear for short times in the upper part of reservoir for 15 years. Nevertheless, interestingly, both approaches yield similar mass transport if the global shape of the solute plume and maximum and minimum values are scrutinized. This comparison of model solutions is better appreciated with breakthrough curves at observation points. As observed in Figure 12, there are differences in the time of

40

ACCEPTED MANUSCRIPT

maximum concentration, and for some periods the upscaled model overestimates the

CE

PT

ED

M

AN US

CR IP T

concentration. However, the general trends of the concentration curves are similar.

AC

Figure 11. Comparison between concentration profiles (mg/L) computed using the Darcy scale governing equations and the solution of upscaled equations using the calculated effective coefficients. Here, time evolution is presented, as well as the maximum and minimum concentration values.

41

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 12. Comparison between breakthrough curves at observation points calculated with

ED

M

the Darcy scale model solution and the upscaled model.

7.1 Comparison of dispersion coefficients with field data

PT

As mentioned in the previous subsection, a strict comparison of effective coefficients with those measured or deduced from field-data is not possible at present for two reasons: i)

CE

field tests addressed to evaluate the solute transport in reservoirs are commonly interpreted using the one-equation model, and furthermore, a unique dispersion tensor is reported in the

AC

literature; ii) in this work, we computed the effective coefficients for layered and lenticular reservoirs; thus, in order to validate them, similar reservoirs must be used for comparison. In practice, the in situ arrangement and lithology of fluid reservoirs can vary significantly according to its geological history, which leads to large heterogeneities of petrophysical, fluid-rock, and solute transport properties. In one attempt to show such heterogeneities, we plot dispersion coefficients estimated from tracer tests in water and petroleum reservoirs in Figure 13. Note the variety of lithologies and the magnitude of coefficient dispersions 42

ACCEPTED MANUSCRIPT

ranging from 1×10-10 to 1×100 m2/s. In our approach, we assumed a continuous facies in the

 -region, and furthermore the coefficient D is dominant. Surprisingly, depending on the rock permeabilities, the numerical estimations of

D 

 xx

can be very close to those

dispersion coefficients estimated for: gravel, layered compost and limestone gravel,

ED

M

AN US

CR IP T

laminated clay, layers of silt and sand, fractured carbonates, and open faults.

PT

Figure 13. Comparison of dispersion coefficients computed in this work [  D  xx ,  D  xx

CE

,  D  xx and  D  xx ], and some reported in the literature (most of the field data are interpreted with a unique dispersion coefficient contained in the advection-dispersion

AC

equation model; there are no data using the two-equation model) for various lithologies in water or petroleum reservoirs. [43]: alternations of sandstone and sand bodies with shale

layers; [44]: layered units: silt, sand, highly weathered and fresh bedrock; [45]: fractured carbonate with tight limestone to vuggy and dolomitized facies; [46]: gravel, including

pebble and boulders; [47]: consecutive layers of compost and a layer of limestone gravel; [48]: karst aquifer; [49]: karst aquifer; [50]: open conductive fault and sand matrix; [51] sand and gravel aquifer with laminated or lensed clay. 43

ACCEPTED MANUSCRIPT

8. Conclusions Applying the volume averaging method, we develop a two-equation reservoir scale model for solute mass transport through two well-defined geological facies. Several effective

CR IP T

coefficients arose, such as: mass dispersion tensors, pseudo convective and pseudo absorption terms, and interfacies mass transfer coefficients. The definitions of these coefficients are provided with the closure problems used to estimate them, which allows us to link them with physics and the geometrical distribution of facies at a smaller scale, i.e., the Darcy scale. For lenticular and layered arrangements of facies, the longitudinal

AN US

components, if applicable, of effective coefficients are computed by solving the associated closure problems as function of rock permeability. On the basis of numerical estimations, we found that mass dispersion tensors and pseudo velocities are the relevant coefficients. The averaged two-equation model is compared with direct simulations at the Darcy scale when a solute slug is injected into a synthetic lenticular reservoir. It was found that general

M

trends of solute plume and breakthrough curves agree well. In this work, we established the fundamentals to the derive a one-equation model, which is widely used to interpret tracer

PT

Acknowledgements

ED

injection tests in hydrocarbon, water, and geothermal reservoirs (ongoing work).

AC

CE

ECHH thanks the support from UANL during his research stay at FTC-Linares.

44

ACCEPTED MANUSCRIPT

References [1] Quintard, M., Whitaker, S. Transport in chemically and mechanically heterogeneous porous media IV: large-scale mass equilibrium for solute transport with adsorption. Advances in Water Resources 1998; 22(1): 33-57.

CR IP T

[2] Porta, G. M., Riva, M., Guadagnini, A. Upscaling solute transport in porous media in the presence of an irreversible bimolecular reaction. Advances in Water Resources 2012; 35: 151-162.

[3] Wood, B. D. The role of scaling laws in upscaling. Advances in Water Resources 2009;

AN US

32(5): 723-736.

[4] Espinosa-Paredes, G. Heat transfer processes upscaling in geoenergy fields. Energy Sources, Part A: Recovery, Utilization, and Environmental Effects 2014; 36(20): 22542262.

[5] Chabanon, M., David, B., Goyeau, B. Averaged model for momentum and dispersion in

M

hierarchical porous media. Physical Review E 2015; 92(2): 023201.

ED

[6] Raoof, A., Hassanizadeh, S. M. Upscaling transport of adsorbing solutes in porous media. Journal of Porous Media 2010; 13(5): 395-408.

PT

[7] Hill, D. J., Minsker, B. S., Valocchi, A. J., Babovic, V., Keijzer, M. Upscaling models of solute transport in porous media through genetic programming. Journal of

CE

Hydroinformatics 2007; 9(4): 251-266. [8] Geiger, S., Cortis, A., Birkholzer, J. T. Upscaling solute transport in naturally fractured

AC

porous media with the continuous time random walk method. Water Resources Research 2010; 46(12): W12530. [9] Sahimi, M. Dispersion in porous media, continuous-time random walks, and percolation. Physical Review E 2012, 85: 016316. [10] Ghanbarian-Alavijeh, B., Skinner, T.E., Hunt, A.G. Saturation dependence of dispersion in porous media. Physical Review E 2012, 85: 066316. 45

ACCEPTED MANUSCRIPT

[11] Marjerrison, D. M., Fassihi, M. R. A procedure for scaling heavy-oil combustion tube results to a field model. In SPE/DOE Enhanced Oil Recovery Symposium. Society of Petroleum Engineers 1992; 22-24 April: Tulsa, Oklahoma. [12] Islam, M. R., Ali, S. F. New scaling criteria for in-situ combustion experiments.

CR IP T

Journal of Petroleum Science and Engineering 1992; 6(4): 367-37. [13] Aguilar-Madera, C.G., Cazarez-Candia, O., Valdés-Parada, F.J. Volume averaged equations for mass transport and reaction for in-situ combustion. International Journal of Chemical Reactor Engineering 2017; 15(5): DOI: 10.1515/ijcre-2017-0124.

[14] Salimi, H., Bruining, J. Upscaling of fractured oil reservoirs using homogenization

AN US

including non-equilibrium capillary pressure and relative permeability. Computational Geosciences 2012; 16(2): 367-389.

[15] Ahmadi, A., Aigueperse, A., Quintard, M. Upscaling of nonwetting phase residual transport in porous media: a network approach. Transport in Porous Media 2001; 43(2):

M

309-353.

[16] Pickering, K., Stow, D., Watson, M., Hiscott, R. Deep-water facies, processes and

ED

models: a review and classification scheme for modern and ancient sediments. EarthScience Reviews 1986; 23(2): 75-174.

PT

[17] Shanmugam, G. 50 years of the turbidite paradigm (1950s—1990s): deep-water processes and facies models—a critical perspective. Marine and Petroleum Geology 2000;

CE

17(2): 285-342.

[18] Cullis, S., Colombera, L., Patacci, M., McCaffrey, W. D. Hierarchical classifications

AC

of the sedimentary architecture of deep-marine depositional systems. Earth-Science Reviews 2018; 179: 38-71. [19] McCarthy, M., Pritchard, H., Willis, I., King, E. Ground-penetrating radar measurements of debris thickness on Lirung Glacier, Nepal. Journal of Glaciology 2017; 63(239): 543-555.

46

ACCEPTED MANUSCRIPT

[20] Mertes, J. R., Thompson, S. S., Booth, A. D., Gulley, J. D., Benn, D. I. A conceptual model of supra‐glacial lake formation on debris‐covered glaciers based on GPR facies analysis. Earth Surface Processes and Landforms 2017; 42(6): 903-914. [21] Watson, C. S., Quincey, D. J., Carrivick, J. L., Smith, M. W., Rowan, A. V., Richardson, R. Heterogeneous water storage and thermal regime of supraglacial ponds on

CR IP T

debris‐covered glaciers. Earth Surface Processes and Landforms 2018; 43(1): 229-241.

[22] Beigi, M., Jafarian, A., Javanbakht, M., Wanas, H. A., Mattern, F., Tabatabaei, A. Facies analysis, diagenesis and sequence stratigraphy of the carbonate-evaporite succession of the Upper Jurassic Surmeh Formation: Impacts on reservoir quality (Salman Oil Field,

AN US

Persian Gulf, Iran). Journal of African Earth Sciences 2017; 129: 179-194.

[23] Marchionda, E., Deschamps, R., Cobianchi, M., Nader, F. H., Di Giulio, A., Morad, D. J., Ceriani, A. Field-scale depositional evolution of the Upper Jurassic Arab Formation (onshore Abu Dhabi, UAE). Marine and Petroleum Geology 2018; 89: 350-369.

M

[24] Li, F., Jing, X., Zou, C., Zhang, H., Xiang, F. Facies analysis of the Callovian– Oxfordian carbonates in the northeastern Amu Darya Basin, southeastern Turkmenistan.

ED

Marine and Petroleum Geology 2017; 88: 359-380. [25] He, Y., Mu, H., Kang, Y., Wang, Y., Fan, B. Characteristics of an Upper Jurassic

PT

Carbonate Ramp in the Northern Amu-Darya Basin. International Journal of Geosciences 2018; 9(2): 148.

CE

[26] Jadoon, Q. K., Roberts, E. M., Henderson, B., Blenkinsop, T. G., Wüst, R. A., Mtelela, C. Lithological and facies analysis of the Roseneath and Murteree shales, Cooper Basin,

AC

Australia. Journal of Natural Gas Science and Engineering 2017; 37: 138-168. [27] Li, Y., Schieber, J., Fan, T., Wei, X. Pore characterization and shale facies analysis of the Ordovician-Silurian transition of northern Guizhou, South China: The controls of shale facies on pore distribution. Marine and Petroleum Geology 2018; 92: 697-718. [28] Song, C., Liu, Z., Wang, Y., Li, X., Hu, G. Multi-waveform classification for seismic facies analysis. Computers & Geosciences 2017; 101: 1-9. 47

ACCEPTED MANUSCRIPT

[29] Zhao, T., Li, F., Marfurt, K. J. Seismic attribute selection for unsupervised seismic facies analysis using user-guided data-adaptive weights. Geophysics 2018; 83(2): O31O44. [30] Strebelle, S. B., Journel, A. G. Reservoir modeling using multiple-point statistics. In

30 September-3 October: New Orleans, Louisiana.

CR IP T

SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers 2001;

[31] Zhi, G., Longde, S., Ailin, J., Tao, L. U. 3D geological modeling for tight sand gas reservoir of braided river facies. Petroleum Exploration and Development 2015; 42(1): 8391.

AN US

[32] Soete, J., Claes, S., Claes, H., Janssens, N., Cnudde, V., Huysmans, M., Swennen, R. Lattice Boltzmann simulations of fluid flow in continental carbonate reservoir rocks and in ppscaled rock models generated with multiple-point geostatistics. Geofluids 2017; ID 7240524, https://doi.org/10.1155/2017/7240524

M

[33] Bayat, A., Asghari, O., Bahroudi, A., Tavvakoli, M. Facies Modeling of Heterogeneous Carbonates Reservoirs by Multiple Point Geostatistics. Journal of

ED

Petroleum Science and Technology 2016; 6(2): 56-65. [34] Hunt, A.G., Skinner, T.E., Ewing, R.P., Ghanbarian-Alavijeh, B. Dispersion of solutes

PT

in porous media. The European Physical Journal B 2011; 80: 411-432. [35] Valdés-Parada, F.J., Lasseux, D., Bellet, F. A new formulation of the dispersion tensor

CE

in homogeneous porous media. Advances in Water Resources 2016; 90: 70-82. [36] Howes, F.A., Whitaker, S. The spatial averaging theorem revisited. Chemical

AC

Engineering Science 1985; 40(8): 1387-1392. [37] Gray, W.G. A derivation of the equations for multiphase transport, Chemical Engineering Science 1975; 30: 229-231. [38] Whitaker S. The method of volume averaging, Kluwer Academic Publishers, The Netherlands 1999.

48

ACCEPTED MANUSCRIPT

[39] Ochoa-Tapia, J.A., Stroeve, P., Whitaker, S. Diffusive transport in two-phase media: Spatially periodic models and Maxwell’s theory for isotropic and anisotropic systems, Chemical Engineering Science 1994; 49: 709-726. [40] Chang, H.C. Multiscale analysis of effective transport in periodic heterogeneous

CR IP T

media, Chemical Engineering Communications 1982; 15: 83-91. [41] NŒtinger, B., Estebenet, T., Quintard, M. Up-scaling flow in fractured media: Equivalence between the large scale averaging theory and the continuous time random walk method, Transport in Porous Media 2001; 43: 581-596.

[42] Valdés-Parada, F.J. Integral formulation for the solution of closure problems in

AN US

upscaling processes, Revista Mexicana de Ingeniería Química 2010; 9: 53-66.

[43] Giadom, F.D., Akpokodje, E.G., Tse, A.C. Determination of migration rates of contaminants in a hydrocarbon-polluted site using non-reactive tracer test in the Niger Delta, Nigeria, Environmental Earth Sciences 2015; 74: 879-888.

M

[44] Jeon-Woo, K., Heechul, C., Jin-Yong, L. Characterization of hydrogeologic properties for a multi-layered alluvial aquifer using hydraulic and tracer tests and electrical resistivity

ED

survey, Environmental Geology 2005; 48: 991-1001. [45] Bogatkov, D., Babadagli, T. Fracture network modeling conditioned to pressure

75: 154-167.

PT

transient and tracer test dynamic data, Journal of Petroleum Science and Engineering 2010;

CE

[46] Maloszewski, P., Herrmann, A., Zuber, A. Interpretation of tracer tests performed in fractured rock of the Lange Bramke basin, Germany, Hydrogeology Journal 1999; 7: 209-

AC

218.

[47] Díaz-Goebes, M., Younger, P.L. A simple analytical model for interpretation of tracer tests in two-domain subsurface flow systems, Mine Water and the Environment 2004; 23: 138-143.

49

ACCEPTED MANUSCRIPT

[48] Goldscheider, N. A new quantitative interpretation of the long-tail and plateau-like breakthrough curves from tracer tests in the artesian karst aquifer of Stuttgart, Germany, Hydrogeology Journal 2008; 26: 1311-1317. [49] Morales, T., Uriarte, J.A., Angulo, B., Olazar, M., Arandes, J.M., Antigüedad, I. Characterization of flow and transport dynamics in karst aquifers by analyzing tracer test

CR IP T

results in conduits and recharge areas (the Egino Massif, Basque Country, Spain): environmental and management implications, Environmental Earth Sciences 2018; 77: 291. [50] Coronado, M., Ramírez-Sabag, J., Valdiviezo-Mijangos, O. Double-porosity model for tracer transport in reservoirs having open conductive geological faults: Determination of

AN US

the fault orientation, Journal of Petroleum Science and Engineering 2011; 78: 65-77.

[51] Yang, Y.S., Lin, X.Y., Elliot, T., Kalin, R.M. A natural-gradient field tracer test for evaluation of pollutant-transport parameters in a porous-medium aquifer, Hydrogeology

AC

CE

PT

ED

M

Journal 2001; 9: 313-320.

50

ACCEPTED MANUSCRIPT

Appendix A Estimation of orders of magnitude We start by estimating the intrinsic average concentration as a function of nonlocal concentrations around the REV centroid x . Such estimation is provided using a Taylor

C

 x y

 C

 x

 y  C



1  y y :  C x 2

Introducing the following estimates 

r  C   O 0  L 

   

AN US

y  C



CR IP T

series expansion as follows:





 r 2  C 0   O  LC1L 



   



(A.1)

x

(A.2)

(A.3)

M

1 y y :  C 2



ED

we can write that

PT

y  C



C

x

1 y y :  C 2

 x



(A.4) x

C



(A.5) x

r0

L,

r02

LC1L

(A.6)

AC

CE

only if the length-scale constraints

are met. In Eq. (A.3), L and LC1 are the characteristic distances and gradients, respectively, for significant changes in concentration. Indeed, such characteristic lengths could be time-dependent, and for the present case, L and LC1 can be the size of the entire reservoir.

51

ACCEPTED MANUSCRIPT

The inequalities (A.6) can be used to get out average quantities from the average operator. For instance, applying the average operator to the spatial decomposition of concentration, we get C  C



 C

(A.7)

C   C



CR IP T

and on the basis of inequalities (A.6), we have  C

(A.8)

and by comparison to definition (10), we have C  0

AN US

(A.9)

To simplify Eq. (28), we perform the following order of magnitude estimates:   D C A ,REV   n C dA   O  L  V  A 

M

  1    D     V  

D C   D C  O  2   L  



ED



AC

CE

PT

  u C



  

u C   O     L 



  u C



D  C





  O  u LC

(A.11)

(A.12)

u C    u C  O      L 



(A.10)

(A.13)

  

(A.14)

D C   O     L 

(A.15)

  

D C  D C  O      L 

(A.16)

Comparing Eqs. (A.10) and (A.11), we realize that 52

ACCEPTED MANUSCRIPT

  1    D     V  

 n C dA     A 



  D C



(A.17)

if the following length-scale constraint is satisfied

  A ,REV  1          V

1

(A.18)

CR IP T

2  L   L

Here, A ,REV is the surface area within the REV. We noticed that this constraint is easier to satisfy if the REV volume is sufficiently larger than the  -region size, i.e.

r0

(A.19)

AN US

L

Here, it must be recalled that L is the characteristic length of the  -region. By inspecting Figure 1, we observe that for layered reservoirs, the inequality A ,REV

1

(A.20)

M

V

ED

is satisfied if a rectangular REV (such as that depicted in the figure) is used.

PT

By comparing the length-scale constraints (A.6) and (A.19), we elucidate that

L

L, LC1

(A.21)

CE

This inequality allows us to conclude that

AC

  u C









  u C ,   u C

D  C







D C

(A.22)

(A.23)

The closure problem is considered as quasi-steady state if the time-scale constraint

D t * 2 L

1

(A.25)

53

ACCEPTED MANUSCRIPT

is satisfied. Here, t * is the characteristic time of the mass transport problem at the Darcy scale. Thus, with inequality (A.25), we can write C t





(A.26)

CR IP T

  D C

Appendix B The closure problems

The closure problems are obtained by substitution of formal solutions of spatial deviations

AN US

[Eqs. (43) and (44)] into governing equations for spatial deviations [Eqs. (38)-(41)]. Thus we develop the following boundary-value problems. B.1. Closure problem 1

M

In the  -region

In the  -region



ED

   u b   u     D b  







PT

  u b    D b 

1 V

1 V



n   u b  D b  dA

(B.1)

A







n  u b  D b dA

(B.2)

A

AC

CE

At the  -boundary

b  b

(B.3)

At the  -boundary





n   u b  D b   n  u b  D b  n  D

(B.4)

B.2. Closure problem 2

54

ACCEPTED MANUSCRIPT

In the  -region









  u b    D b 

1 V







n  u b  D b dA

(B.5)

A

In the  -region







1 V



 A

At the  -boundary

b  b









n  u b  D b  n  u b  D b  n  D

(B.8)

M

B.3. Closure problem 3

(B.6)

(B.7)

AN US

At the  -boundary



n  u b  D b dA

CR IP T



  u b  u    D b 

ED

In the  -region

   u d     u     D d  

1 V

 n   u d  D d  dA A

PT

1  V

(B.9)

 n  u dA A

CE

In the  -region









AC

  u d    D d 

1 V

 n  u d  D d  dA

(B.10)

A

At the  -boundary d  d

(B.11)

At the  -boundary

55

ACCEPTED MANUSCRIPT





n   u d  D d   n  u d  D d  n  u

(B.12)

B.4. Closure problem 4









  u d    D d 

1 V





In the  -region











n  u d  D d dA

A

A

1  V

(B.14)

 n  u dA A

M

At the  -boundary



ED

d  d

At the  -boundary





(B.13)

 n  u d  D d  dA

1 V

AN US

  u d    u    D d 

CR IP T

In the  -region

(B.15)



(B.16)

CE

PT

n  u d  D d  n  u d  D d  n  u

B.5. Closure problem 5

AC

In the  -region

   u s      D s  

1 V



n   u s  D s  dA

(B.17)

A

In the  -region

56

ACCEPTED MANUSCRIPT









  u s    D s 

1 V







n  u s  D s dA

(B.18)

A

At the  -boundary

s  1  s At the  -boundary



CR IP T

(B.19)

n   u s  D s   n  u s  D s



(B.20)

These closure problems must be solved once and uncoupled from the closed volume averaged equations for each facies. With the solution of closure variables, all effective

AN US

coefficients can be calculated and used in the average equations. Nevertheless, the solution of closure problems can be a demanding task in some cases. For this problem, each closure problem consists of a steady-state boundary value problem with two coupled dependent variables; note that the dependent variables are coupled through the inter-facies boundary conditions. For closure problems 1-4 presented above, the eight dependent variables are

M

vectors, whereas closure problem 5 has two scalar dependent variables. For instance,

ED

closure problem 1 contains the closure variables b and b , whose individual components for 2D Cartesian coordinates are

 b  x ,  b  y ,

 b  x

and

b  y .

PT

Furthermore, the efforts to find their solutions become overwhelming as the dimensionality of the problem increases. Additionally, the shape of each partial differential equation lies

CE

within the advection-diffusion-reaction form (see for instance Eq. (B.1)] but containing one integral source which confers nonlocal characteristics.

AC

Within the volume averaging framework [38], it is usual that the closure variables are subjected to periodic conditions at the external cell boundaries, and the problem is wellposed by adding the restriction to the average:

  0,   any closure variable which results from using Eq. (23). However, the usage of homogeneous Dirichlet boundary conditions has demonstrated to yield similar results to those using periodic conditions [39], 57

ACCEPTED MANUSCRIPT

and with the advantage of decreasing the computational or analytical costs. In our calculations, homogeneous Dirichlet boundary conditions were used at the external boundaries of cells. Finally, we mention that the usage of the finite element method to solve the closure problems is not a unique manner of solution. The solution of these problems have been addressed using other techniques, for instance: analytical techniques have been

CR IP T

employed considering regular cells [40], continuous time random walk methods in fractured media [41], and pseudo-analytic procedures using Green functions in reactive mass transport [42].

AN US

Appendix C

Obtaining the closed volume-averaged equation for  -facies

In this section, we exemplify the algebraic procedure to obtain the closed volume average equation for the  -facies. To this end, let us substitute the formal solution (43) into the

AC

CE

PT

ED

M

average equation (25) to get

58

ACCEPTED MANUSCRIPT



1 V



     u

t

 n   u

C





    D    C   C  dA

C

 D















A

  1    D    V  

 n b dA      C A 

  1    D    V     1    D    V  

 n d dA     C A 



  1     D     V   





 n b dA      C A 

   1     D    n d dA  C   V   A           n s dA   C  C  A       u b  C     u b  C     u d C           u d C      u s  C  C          1    n   u b  D b  dA  C  V   A     1   n  u b  D b dA  C  V   A    1     n   u d  D d  dA  C V   A     1    n  u d  D d dA C V   A     1    n   u s  D s  dA  C  C   V    A 

AN US





   







   



M





CR IP T





 C

PT

ED





AC

CE





(C.1)

Based on the length-scale constraints presented in the Appendix A, the quantities  C  C



, C



, C



and  C 



 C





,

 are treated as constants within the integrals  

59

ACCEPTED MANUSCRIPT

of area

   dA and volume  

, and furthermore, the terms can be simplified. The next

A

step involves the grouping of similar terms in Eq. (C.1) to get

t





     u

 1      D  D   V 

 1    D   V 

CE AC



  n b dA  u b    C     A 

 n b dA 

u b

A

 n d dA  A

 n d dA  A

 n s dA  A

   C      C   



u d

u d

 





     C      

   C   

u s





   



 C



    

ED

   n   u d  D d  u  dA  C A    n  u d  D  d dA        C A    n  u b  D  b  D dA            C A     C  n  u b  D  b dA        A       n   u s  D s  dA  C  C  A 

PT

 1  V   1  V   1  V   1  V   1  V 



AN US

 1    D   V 

C

M

 1    D   V   1    D   V 



CR IP T



 C









(C.2)

60

ACCEPTED MANUSCRIPT

From this equation, is already possible to define the several effective coefficients presented in the closed average  -facies equation (46). A similar theoretical procedure can be outlined for the equation for the  -facies [Eq. (47)], but such algebraic effort is omitted for

AC

CE

PT

ED

M

AN US

CR IP T

the sake of brevity.

61

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

Graphical abstract

62