PoreFlow: A complex pore-network model for simulation of reactive transport in variably saturated porous media

Author's Accepted Manuscript

PoreFlow: A complex pore-network model for simulation of reactive transport in variably saturated porous media A. Raoof, H.M. Nick, S.M. Hassanizadeh, C.J. Spiers

www.elsevier.com/locate/cageo

PII: DOI: Reference:

S0098-3004(13)00226-4 http://dx.doi.org/10.1016/j.cageo.2013.08.005 CAGEO3249

To appear in:

Computers & Geosciences

Received date: 13 December 2012 Revised date: 3 July 2013 Accepted date: 19 August 2013 Cite this article as: A. Raoof, H.M. Nick, S.M. Hassanizadeh, C.J. Spiers, PoreFlow: A complex pore-network model for simulation of reactive transport in variably saturated porous media, Computers & Geosciences, http://dx.doi.org/ 10.1016/j.cageo.2013.08.005 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 galley proof before it is published in its final citable 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.

1 2 3

PoreFlow: A Complex Pore-Network Model for Simulation of Reactive Transport in Variably Saturated Porous Media A. Raoof

4

7

, H.M. Nicka,b , S.M. Hassanizadeha , C.J. Spiersa∗

a

5 6

a∗

b

Department of Earth Sciences, Utrecht University, Utrecht, The Netherlands Department of Geoscience and Engineering, Delft University of Technology, Delft, The Netherlands

8

Abstract

9

This study introduces PoreFlow, a pore-network model capable of simulat-

10

ing fluid flow and multi-component reactive and adsorptive transport under

11

saturated and variably saturated conditions. PoreFlow includes a variety of

12

modules, such as: pore network generator, drainage simulator, calculation of

13

pressure and velocity distributions, and modeling of reactive solute transport

14

accounting for advection and diffusion. The pore space is represented us-

15

ing a multi-directional pore-network capable of capturing the random struc-

16

ture of a given porous media with user-defined directional connectivities for

17

anisotropic pore structures. The reactions can occur within the liquid phase,

18

as well as between the liquid and solid phases which may result in an evolu-

19

tion of porosity and permeability. Under variably saturated conditions the

20

area of interfaces changes with degree of the fluid saturation.

21

PoreFlow uses complex formulations for more accurate modeling of transport

22

problems in presence of the nonwetting phase. This is done by refining the

23

discretization within drained pores. An implicit numerical scheme is used ∗

Corresponding author E-mail address: [email protected]

Preprint submitted to Computers and Geosciences

September 4, 2013

1

to solve the governing equations, and an efficient substitution method is ap-

2

plied to considerably minimize computational times. Several examples are

3

provided, under saturated and variably saturated conditions, to demonstrate

4

the model applicability in hydrogeology problems and petroleum fields. We

5

show that PoreFlow is a powerful tool for upscaling of flow and transport in

6

porous media, utilizing different pore scale information such as various inter-

7

faces, phase distributions and local fluxes and concentrations to determine

8

macro scale properties such as average saturation,

9

solute dispersivity, adsorption coefficients, effective diffusion and tortuosity.

10

Such information can be used as constitutive relations within continuum

11

scale governing equations to model physical and chemical processes more

12

accurately at the larger scales.

13

Keywords: pore scale, pore network modeling, variably saturated, reactive

14

transport

15

1. Introduction

c1

relative permeability,

Simulation of variably-saturated flow and (multi-component) reactive trans-

16 17

port is highly important in many applications involving both natural (soils,

18

rocks, woods, hard and soft tissues, etc.) and industrial (concrete, bioengi-

19

neered tissues, etc.) porous media. Compared to studies of solute transport

20

at the macro scale, e.g., column scale or field scale (Class et al., 2008; Nick

21

et al., 2011; Geiger et al., 2012), there is less comprehensive research on this

22

subject at the pore scale (Raoof and Hassanizadeh, 2009; Bijeljic et al., 2011;

23

Ovaysi and Piri, 2011). Study of pore-scale processes is essential for explana-

24

tion and understanding of flow and transport processes at the macro-scale. In c1

retaliative

2

1

addition, pore scale experimental works are often difficult and/or expensive

2

to perform; this makes the use of pore-scale modeling very appealing.

3

1.1. Issue of scales

4

Transport in porous media is a multi-scale problem with each scale con-

5

taining specific information about the underlying physical process (Kechagia

6

et al., 2002; Held and Celia, 2001). The scale hierarchy associated with flow

7

and transport problems in porous media is often divided into: molecular

8

scale, micro or pore scale, macro or lab scale, meso or field scale, and mega

9

or regional scale.

10

While modeling at a larger scale, it is usually not feasible to take all pore-

11

scale properties, such as interfaces, into account. However, their lumped ef-

12

fects must be included in macro-scale descriptions (Raoof and Hassanizadeh,

13

2010).

14

macro scale representation of processes, it is a useful tool to study averaging

15

effects, to explore how pore scale processes manifest themselves in macro-

16

scopic level (core or field scale). Pore-scale models can be used to elucidate

17

the pore-scale nature of processes such as: spreading of solute due to deviac1

c1

Since pore-scale modeling provides a bridge between pore scale and

Pore-scale models can be used to elucidate the pore-scale nature of processes such as:

spreading of solute due to deviations in pore velocities, adsorption of solutes or colloids to the grain surfaces under saturated conditions, reaction at the nonwetting-wetting interface under variably saturated conditions. Since pore-scale modeling provides a bridge between pore scale and macro scale representation of processes, it is a useful tool to study averaging effects. While this is true, we do not use pore scale models for calculating average concentration that is a step in finding relationships between macro-scale and pore-scale parameters and properties.

3

1

tions in pore velocities, adsorption of solutes or colloids to the grain surfaces

2

under saturated conditions, reaction at the nonwetting-wetting interface un-

3

der variably saturated conditions.

4

1.2. Pore-network modeling applications

5

Modeling fluid flow in porous media at the macro scale requires specifi-

6

cation of the constitutive properties of the porous medium. Examples are:

7

i) the relationship between the capillary pressure, Pc (the difference between

8

the pressures of the nonwetting and wetting fluids) and the fluid saturation,

9

Sw , ii) the relative permeability, kr , as a function of either the saturation

10

or capillary pressure, and iii) the relationship between solute dispersivity or

11

effective diffusion and saturation in solute transport processes.

12

The relative permeability of a fluid is a measure of the conductance of the

13

porous medium for that fluid at a given saturation. Relative permeability

14

measurements on field samples are difficult and time consuming. In gen-

15

eral, experimental determination of the Pc − Sw relationship is easier than

16

measurement of relative permeability. For this reason, empirical relation-

17

ships, such as those of Brooks and Corey (Brooks and Corey, 1964) and Van

18

Genuchten (Van Genuchten, 1980), are often used to model the dependence

19

of relative permeability on capillary pressure or saturation (e.g. Nick and

20

Matth¨ai, 2011).

21 22

c1

As another approach, one may use pore network models (PNMs) to obtain

constitutive properties. One of the early attempts to estimate relative permec1

Another approach for obtaining constitutive properties is to use Pore-Network Models

(PNMs).

4

1

ability was to use bundle-of-capillary-tubes models. This was based on the

2

assumption that a porous medium may be modeled as bundle of capillary

3

tubes of various diameters. However, such models ignore the interconnected

4

nature of porous media and often do not provide realistic results. PNMs,

5

first suggested by (Fatt, 1956), offer a more realistic approach for calculating

6

constitutive properties. The vast majority of PNMs consist of pore bodies (or

7

nodes) and pore throats (or channels), along with a specified topological con-

8

figuration which prescribes how pore bodies are connected via pore throats.

9

Pore bodies are meant to represent larger void spaces found in natural porous

10

media. The narrow openings that connect adjacent pore bodies are modeled

11

by pore throats, which are essentially capillary tubes. The pore-network

12

approach for modeling multiphase flow properties has been employed exten-

13

sively in the petroleum engineering literature (Chatzis and Dullien, 1985;

14

Celia et al., 1995; Blunt, 2001; Algive et al., 2012)c2 . In recent years, various

15

pore network modelings have been developed to explore chemical and bio-

16

logical processes, such as: dissolution of organic liquids (Zhou et al., 2000;

17

Held and Celia, 2001; Knutson et al., 2001); CO2 sequestration (Kang et al.,

18

2010; Kim et al., 2011; Mehmani et al., 2012; Raoof et al., 2012; Varloteaux

19

et al., 2012)c3; biomass growth (Suchomel et al., 1998a; Kim and Fogler,

20

2000; Dupin et al., 2001; Gharasoo et al., 2011; Rosenzweig et al., 2013);

21

solute dispersivity (Vasilyev et al., 2012)c4 , and adsorption (Acharya et al.,

22

2005; Li et al., 2006; Raoof and Hassanizadeh, 2010; K¨ohne et al., 2011)c5 . c2 c3 c4 c5

New New New New

reference reference reference reference

added added added added

5

1

Because of their ability to simulate the highly disordered geometry of pore

2

space and relatively low computational cost (e.g. in comparison with Lat-

3

tice Boltzmann simulation, Sholokhova et al., 2009)c1 , PNMs hold promise

4

as tools for predicting multiphase flow properties of specific porous media.

5

For example, the dependence of capillary pressure on saturation is modeled

6

by determining the location of fluid-fluid interfaces throughout the network

7

using the Young-Laplace equation. This is sometimes modified by other pore-

8

level mechanisms, such as snapoff during imbibition (e.g., Yu and Wardlaw,

9

1986). Also, the dependence of relative permeability on saturation is de-

10

termined by computing the resistance to flow in the connected portion of

11

a fluid. In these calculations, resistance to the flow within the pore bodies

12

is commonly ignored, assuming that conductance within the pore bodies is

13

much larger than the conductance within the pore throats. Commonly, an

14

average pressure and/or concentration is assigned to a given pore body or

15

pore throat. For each pore, change in solute mass is described by mass bal-

16

ance equations (Lichtner, 1985). Using information on local surface area, and

17

applying an area-normalized reaction rate, a kinetic reaction is calculated for

18

each pore (Raoof and Hassanizadeh, 2010).

19

One shortcoming of a pore-network model is its idealization of the pore

20

space using simple geometries; often, pores are assumed to have uniform

21

circular or square cross-sectional shapes. This makes it difficult to simu-

22

late some processes, such as biogeochemical reactions, that could lead to

23

a significant change in complex pore geometries as a result of dissolution, c1

New reference and text added

6

1

precipitation, and/or biological clogging. However, pore-network models

2

are computationally cheap, and recent advances have allowed modeling a

3

degree of irregularity in channel cross-sectional shapes. In addition, pore-

4

network models are capable of capturing important statistical characteristics

5

of porous media such as pore size distributions (Oren et al., 1998; Lindquist

6

et al., 2000), together with coordination number distributions (Raoof and

7

Hassanizadeh, 2009) and topological parameters such as Euler number (Vo-

8

gel and Roth, 2001).

9

1.3. Objective

10

The objective of this article is to introduce a simulation tool, called Pore-

11

Flow, for simulation of pore-scale flow and reactive/adsorptive transport in

12

two- and three-dimensional domains. PoreFlow has several features which

13

makes it a powerful tool for pore-scale modeling. It allows us to model sev-

14

eral different processes, starting from generation of a pore network model for

15

a specific porous media up to the simulation of

16

fluid flow and c2 reactive/adsorptive solute or colloid transport (with possible

17

changes in pore sizes), under both saturated and variably saturated condi-

18

tions.

19

throats and, unlike other pore network simulators, it refines discretizations of

20

drained pores to simulate flow and transport more accurately in the presence c1 c2 c3

c3

c1

drainage and (two-phase)

In addition, PoreFlow assigns volumes to both pore bodies and pore

Text added. ractive In addition, unlike other pore network simulators, PoreFlow refines discretizations

of drained pores to simulate flow and transport more accurately in the presence of the nonwetting phase.

7

1

of the nonwetting phase. Doing so, we can calculate velocity and concentra-

2

tion distributions within drained pores. Using such a model, one can obtain

3

different pore-scale information for a given porous structure and derive vari-

4

ous constitutive relations to be used within macro-scale models appropriate

5

for modeling at the continuum scale.

6

We start with a brief description of the software structure, followed by

7

an explanation of the model features. This includes a presentation of the

8

governing equations, and their discretization in space and time. Then, several

9

examples are used to illustrate applicability of this code under various flow

10

and transport conditions.

11

2. PoreFlow structure and methods In this section, we present the structure of the modeling tool, simulation

12 13

modes, and the coupling between different processes.

14

2.1. Coupled model

15

PoreFlow uses a Multi-Directional Pore-Network (MDPN) to represent

16

porous media, and performs flow and solute transport simulations in such a

17

network. The numerical scheme is developed based on the assumption that

18

porous medium can be represented as a network of pore bodies and pore

19

throats, both with finite volumes. There are two options to simulate the

20

transport part: i) advection only, ii) advection and diffusion. When includ-

21

ing multicomponent chemical reactions, Biogeochemical Reaction Network

22

Simulator (BRNS) (Regnier et al., 2002)c1 performs the adsorption/reaction c1

reference added.

8

1

part. This gives major advantages for simulation of complex reaction sys-

2

tems.

3

In the case of chemical reactions with the solid surface (for example, dis-

4

solution/precipitation processes), the pore sizes will change, which in turn

5

will cause changes in porosity and permeability. In such cases, after each

6

reaction time step the flow field is updated. In the discretized time domain,

7

the solute transport and reaction processes are solved using a non-iterative

8

sequential split operator approach. The model benefits from an adaptive

9

time stepping approach for geochemical computations within the iterative

10

schemes.

11 12

A summary of computational features of PoreFlow is as follows:

13

– The topology of the pore space is modeled using a Multi-Directional

14

Pore Network (MDPN) which allows a distribution of coordination

15

numbers ranging between one and 26.

16

– To take into account the angularity of pores in natural porous media,

17

pore throats are assigned a variety of cross sectional shapes, with a

18

wide range of shape factors. This includes rectangular, circular, and

19

various irregular triangular cross sections.

20

– Both pore bodies and pore throats have volumes. This means we solve

21

mass balance equations and calculate solute concentrations and mass

22

fluxes within both pore bodies and pore throats.

23

– The pore body sizes can be taken from various statistical distributions.

24

The pore-throat sizes can be either correlated to the pore body size 9

1

distribution or chosen to be independent from pore body sizes.

2

– A drained pore body is refined into smaller compartments occupied

3

by the wetting phase, each with its own concentration. This takes

4

into account the effect of limited mixing due to the partial filling of a

5

drained pore body by the nonwetting phase.

6

– Upon invasion of an angular pore throat by the nonwetting phase, each

7

edge of the pore throat will be considered as a separate domain with

8

its own flow rate and concentration.

9

– Fully implicit numerical schemes are employed for calculating solute

10

concentrations at each time step. Efficient substitution method is ap-

11

plied which considerably reduces the computational time.

12 13

– BRNS is used for simulating various multicomponent geochemical and microbiological processes.

14

– Various parameters and relations, such as coefficient of variation of

15

pore velocities, capillary pressure-saturation (P c − Sw ) curves, rela-

16

tive permeability-saturation curves, interfacial area, solute dispersivity-

17

saturation relation, effective diffusion, tortuosity, and the fraction of

18

percolating saturated pores are also computed.

19

Fig. 1 shows the flowchart of PoreFlow. Besides the capability of Pore-

20

Flow to generate Multi-Directional Pore Networks, it is possible to generate

21

pore networks with regular patterns as well as user-defined pore connections

22

(Raoof, 2011). The structure of the network can be optimized (using a Ge-

23

netic Algorithm routine (Raoof and Hassanizadeh, 2009)) in order to obtain 10

1

a pore network with a specific coordination number distribution. During the

2

simulation, the data are transferred to MATLAB, using MATLAB Engine

3

functions, for the purpose of analysis and postprocessing.

4

2.2. Pore-network modeling In this section we introduce properties of the pore network following by

5 6

discretization method used within PoreFlow.

7

2.2.1. Coordination number distribution in MDPN

8

A pore body coordination number, z, is defined as the number of pore

9

throats connected to it. Most pore-network models have a regular network

10

with a fixed coordination number of six for all pore bodies (except those at

11

the boundaries). This means that any given pore body is connected to six

12

neighboring nodes via tubes that are oriented along the lattice axes in three

13

principal directions (directions number 1, 2, and 3 in Fig. 2). However, there

14

is overwhelming evidence that a very wide range of coordination numbers ex-

15

ists in a real porous medium (Jivkov et al., 2013)c1 . For example, Ioannidis

16

et al. (1997) determined the mean coordination number, z, from serial sec-

17

tions of a sandstone core, and found it to be 3.5. Oren and coworkers (Oren

18

et al., 1998; Oren and Bakke, 2003) developed a process-based reconstruction

19

procedure, which incorporated grain size distribution and other petrograph-

20

ical data obtained from 2D thin sections. They reported mean coordination

21

numbers of 3.5 − 4.5 (Oren and Bakke, 2003).

22

c1

Reference added

11

1

One of the main features of MDPN is that pore throats can be oriented

2

in 13 different directions, allowing a maximum coordination number of 26,

3

as shown in Fig. 2. Then, to get a desired coordination number distribution,

4

we follow an elimination procedure to rule out some of the connections. The

5

elimination procedure is such that a pre-specified mean coordination number

6

can be obtained. A coordination number of zero means that the pore body is

7

eliminated from the network, so there is no pore body located at that lattice

8

point. Pore bodies with the coordination number of one could be eliminated

9

when dead end pores are not considered, except if they are located at the

10

inlet or outlet boundaries (so they belong to the flowing fluid backbone). This

11

network generation method was successfully employed to develop networks

12

that very closely modeled real sandstone as well as granular media, using

13

information on their coordination number distributions. The detail of the

14

network generation method can be found in (Raoof and Hassanizadeh, 2009).

15

The distribution of coordination number and a representative domain of the

16

network used in this study are given in Fig. (3).

17

FIGURE

18

FIGURE

19

2.2.2. Determination of the pore cross section and corner half angles

20

A key characteristic of real porous media is the angular form of pores. It

21

has been demonstrated that having pores with a circular cross section, and

22

thus single-phase occupancy, causes insufficient connectivity of the wetting

23

phase and as a result poor representation of experimental data (Zhou et al.,

24

2000). Angular cross sections retain the wetting fluid in their corner and allow

25

two or more fluids to flow simultaneously through the same pore. Pores which 12

1

are angular in cross section are thus a much more realistic model than the

2

commonly employed cylindrical shape. In PoreFlow, different pore throats

3

are assigned a variety of cross sectional shapes including circular, rectangular,

4

and scalene triangular. The shape of an angular pore cross section is prescribed in terms of a dimensionless shape factor, G (Mason and Morrow, 1991), which is defined as: G=

A , P2

(1)

5

where A and P are the area and the perimeter of the cross section, respec-

6

tively. The shape factor replaces the irregular and complicated shape of a

7

pore throat by an equivalent either circle, triangle, or rectangle shape. The

8

value of shape factors range from a minimum of zero corresponding to a

9

slit to a maximum of 0.08 corresponding to a circular cross section. Values

10

of shape factors for triangular cross sections vary from zero to 0.048 (with

11

maximum for equilateral triangle), and rectangular cross section from zero

12

to 0.062 (with maximum for square shape).

13

2.2.3. Pore space discretization

14

PoreFlow uses a new approach to c1 discretizing the pore space for calcula-

15

tion of the fluid flows and solute concentrations. The continuous pore space

16

is

17

tion, and adsorbed mass, etc. The discretization depends on saturation state

18

of a pore as follows.

c2

c1 c2

discretized into “pore elements”, each with its own pressure, concentra-

descretizing descretized

13

1

In the case of a saturated pore, one pressure and one concentration values

2

are assigned for each pore element representing either a saturated pore body

3

or a pore throat. In the presence of the nonwetting phase, a given pore body

4

or pore throat could be invaded and partially filled by the nonwetting phase,

5

forcing the wetting phase to flow only along the edges. Fig. 4a shows a

6

schematic example of two pore bodies connected to each other using a pore

7

throat with an angular cross section. Under such conditions, pore throats

8

are not necessarily the narrowest constriction along the flow path, anymore.

9

In fact, resistance to the flow within the edges of drained pore bodies may

10

be comparable to, or even larger than, the resistance to the flow within the

11

pore throats. The significance of these effects on flow and solute transport

12

has been addressed by Raoof and Hassanizadeh (2011, 2013)

13

culating, and making comparison with the experimental results, the relative

14

permeabilities and solute dispersivities by considering conductances, and cal-

15

culating different fluxes, within the pore bodies as well as pore throats. It is

16

worth mentioning that, the lower conductivity of a drained pore body also

17

reduces the connectivity and therefore the solute mass flux among saturated

18

pore throats connected to a drained pore body (such as pore throats number

19

2, 3, 4 in Fig. 4b). Under such conditions, a drained pore body is divided

20

into corner units, each of them being a pore element. As shown in Fig 4a,

21

a corner unit is composed of a corner together with half pore body edges

22

connected to it. For example, a cubic pore body consists of 8 different pore c3

c3

through cal-

through calculating the relative permeability-saturation relationship by considering

conductances, and calculating different fluxes, within the pore bodies as well as pore throats.

14

1

elements with 8 different pressure and concentration values assigned to them.

2

Fluid flow and solute mass fluxes between these elements occurs through 12

3

edges within drained pore body.

4

In a drained pore throat with angular cross section, there will be flow

5

of wetting phase only through the edges of the pore throats. Since in an

6

angular pore edges may have different angles, they have different flow rates,

7

and presence of the nonwetting phase limits mixing between them. To take

8

this effect into account, each edge of a drained pore throat is treated as a

9

pore element with its own flow rate and concentration value.

10

The conductance of each edge within drained pore bodies or pore throats

11

needs to be determined as a function of the thickness of the wetting film

12

residing in the edge. This thickness depends on the radius of curvature of

13

the fluid-fluid interface, which in turn depends on the capillary pressure.

14

Corner units of a given drained pore body are connected to the neighboring

15

pore body corners via pore throats. This means that we need to specify

16

connections of pore throats to the corners of pore bodies, which is done

17

through a random process. The detail of algorithm which has been used

18

to associate different pore throats to different corners of neighboring pore

19

bodies can be found in (Raoof and Hassanizadeh, 2011).

20

2.3. Variably saturated flow modeling

21

To simulate a drainage process, the nonwetting phase is considered to

22

enter the network through an external reservoir which is connected to the

23

inlet side of the network. The displaced phase escapes through the outlet

24

face on the opposite side.

15

1

2.3.1. Drainage simulation

2

At low flow rates, the progress of the displacement is controlled by cap-

3

illary forces. This forms the basis for the invasion-percolation algorithm

4

(Wilkinson and Willemsen, 1983; Chandler et al., 1982) used to model drainage

5

in PoreFlow. At every stage of the process, nonwetting phase invades all ac-

6

cessible pore bodies and throats with the lowest entry capillary pressure. The

7

entry capillary pressure is given by Laplace equation (Bear, 1988),  Pc = Pn − Pw = γwn

1 1 + r1 r2

 =

2γwn , r∗

(2)

8

where Pw and Pn are the pressure of the wetting and the nonwetting phases,

9

respectively, r ∗ is the mean radius of curvature and γwn is the surface tension.

10

For a capillary tube of radius r, we have r ∗ = r/ cos(θ) (Young-Laplace’s

11

equation), in which θ is the contact angle between fluid interface and capillary

12

wall. The invading fluid enters and fills an available pore throat only when

13

the injection pressure is equal to or larger than the entry capillary pressure

14

of the pore.

15

The wetting phase is hydraulically connected through filaments within edges

16

of angular pores. The capillary pressure is increased incrementally so that

17

fluid-fluid interfaces will move only a short distance before coming to rest in

18

equilibrium at the opening of slightly smaller pore throats.

19

2.3.2. Fluid flow within drained pores

20

After performing the drainage step, knowing the curvature of interfaces

21

and local saturations in each pore, PoreFlow calculates conductances of edges

22

of the drained pores throats and drained pore bodies as well as conductances

16

1

of saturated pores. Conductances are needed to calculate the flow across the

2

network, which is due to the flow of wetting c1 phase in saturated pores as well

3

as along edges of drained pore bodies and pore throats. The conductance

4

of an angular drained pore depends on its degree of local saturation, which

5

is directly related to the radius of curvature of the meniscus formed along

6

the pore edges. Raoof and Hassanizadeh (2011) used a numerical solutions

7

to calculate the dimensionless conductances of drained pores with angular

8

cross section. This was done by numerically solving the dimensionless form

9

of the Navier-Stokes equations and the equation of conservation of mass.

10

Through this method, knowing the geometry of the pore, one can calculate

11

the conductance of a given edge of a pore body or pore throat.

12

3. Simulating flow and transport within the network

13

3.1. Flow simulation The flow field is established in the network by imposing a

14

c1

pressure

15

difference across the network. All other boundaries of the network parallel

16

to the overall flow direction are treated as no-flow boundaries. We assume

17

that the discharge through a given drained pore throat can be prescribed by

18

Hagen-Poiseuille equation: Nedge

qij,tot =



gij,k

k=1

c1 c1

pahse pressures

17

(pj − pi ) , lij

(3)

1

where qij,tot is the total volumetric flow rate through pore throat ij, gij,k is

2

conductance of k th edge of angular pore throat ij, and pi and pj are pressures

3

at pore elements i and j, respectively. Eq. (3) is valid for laminar flow in

4

a wide range of Reynolds numbers and is assumed to be appropriate for

5

describing flow in pores (Bear, 1972). For incompressible, steady-state flow, the sum of discharges into and out

6 7

c1

of a pore body, or a pore-body corner unit in the case of a drained pore

8

body, must be zero. For the pore bodies, the continuity equation may be

9

written as: CU,i Nedge

 n=1

qi,n +

zi 

qij,tot = 0;

j = 1, 2, ..., zi ,

(4)

j=1

10

CU,i shows the number of edges through which corner unit i, within where Nedges

11

a drained pore body, is connected to other corner units, n, within the same

12

pore body. zi is the coordination number of pore i. Eq. (4) is applied to

13

all pore bodies except those on the two flow boundaries where pressures are

14

specified. Combination of Eq. (3) and (4) for all pores results in a linear

15

system of equations, with a sparse, symmetric and positive-definite coefficient

16

matrix, to be solved for pore body pressures. Then, the flow velocity in all

17

pore throats can be calculated using Eq. (3). Considering the network as an

18

REV, the average pore velocity v can be determined as: v=

c1

Qtot L , Vf

Text added.

18

(5)

1

where Qtot is the total discharge through the network which is the sum of

2

fluxes through all pore throats at the inlet or outlet boundary of the network,

3

L is the network length, and Vf is total fluid volume.

4

Knowing the total discharge through the network, we can also calculate the

5

relative permeability. The relative permeability of the network to the wetting

6

phase at a given saturation and capillary pressure is calculated using Darcy’s

7

law, krw =

μw Qtot , k A ΔP/L

(6)

8

where μ is viscosity, k is intrinsic permeability, and ΔP is the pressure dif-

9

ference between inflow and outflow reservoirs. Repetition of this process at

10

consecutively larger imposed capillary pressures results in a graph of capillary

11

pressure versus saturation and relative permeability versus saturation.

12

3.2. Adsorptive solute/colloids transport, without diffusion

13

As mentioned earlier, we assume that the solutes or colloids maybe ad-

14

sorbed to the fluid-fluid (i.e., nonwetting-wetting) interfaces (shown by NW)

15

in drained pores as well as to solid-fluid (i.e., solid-wetting) interfaces (shown

16

by SW). Based on the definition of pore element given in Section 2.2.3, the

17

unknowns to be solved for are either concentrations of saturated pore bodies,

18

ci , and saturated pore throats, cij , or concentrations of edges of drained pore

19

throats, cij,k , and corner units of drained pore bodies, cCU,i . The adsorbed

20

mass concentrations in the case of saturated pores will be the adsorbed mass

21

concentrations of saturated pore bodies, ssw i , and saturated pore throats,

22

ssw ij . In the case of variably saturated pores, they will be adsorbed mass 19

1

concentrations at SW and NW interfaces of drained pore throats, ssw ij,k , and,

2

snw ij,k , respectively, and adsorbed to interfaces of corner units of drained pore

3

nw bodies, ssw CU,i , and, sCU,i , respectively.

4

To present the formulation, we provide the mass balance equations for a

5

pore body corner unit and a drained angular pore throat, similar to those

6

shown in Fig. 4a. Indeed, the mass balance equations should be written for a

7

given chemical component, l, with the concentration of cl . However, to keep

8

the equations less crowded we drop the superscript l from the concentrations

9

throughout this manuscript.We assume that the flow is from corner unit j

10

towards corner unit i through corners of drained pore throat ij. For a given

11

corner unit (with concentration cCU,i and volume VCU,i ), we can write the

12

following mass balance equation: throat

N ij

N CU,i

Nin edge in,edge    d VCU,i (cCU,i ) = qij,k cij,k + qi,n cn − QCU,i cCU,i dt n=1 j=1 k=1

−VCU,i

d  sw  d  nw  sCU,i − VCU,i s , dt dt CU,i

(7)

13

ij where the first term on the r.h.s. is due to the mass influx via Nedge edges

14

throat of Nin throats with flow towards the corner unit i. The second term on

15

CU,i the r.h.s. accounts for the mass influx from Nin,edge neighboring corner units

16

(within the same pore body) with flow towards corner unit i. The third term

17

shows the mass leaving the corner unit, where QCU,i is the total flux leaving

18

(or entering) the corner unit i. The last two terms account for the mass

19

adsorbed onto the solid walls of the corner unit and on the NW interface

20

−3 nw within the corner unit, respectively. ssw CU,i [ML ] and sCU,i are the mass

21

adsorbed to SW and NW interfaces per unit volume of the corner unit. 20

1

We should note that, in the case of saturated pores, the second term on the

2

ij right-hand-side of Eq. (7) vanishes, and the value of Nedge in the first term

3

will be equal to one, since there is no edge flow present.

4

nw To close the system, we need extra equations for ssw CU,i and sCU,i . Here, we

5

assume linear equilibrium adsorption at both SW and NW interfaces, such

6

that:

sw sw ssw CU,i = kd,i aCU,i cCU,i nw nw snw CU,i = kd,i aCU,i cCU,i ,

(8)

sw nw where kd,i and kd,i [L] are pore scale adsorption distribution coefficients at

SW and NW interfaces, respectively. The specific surface area, aαw CU,i , is defined as

Aαw CU,i = , (9) VCU,i is the total area of the αw interface within corner unit i (where aαw CU,i

7

where Aαw CU,i

8

α = s, n).

9

Next, the mass balance equation for the k th edge element of a drained

10

pore throat may be written as (assuming that corner unit j is the upstream

11

node):

Vij,k

dcij,k d  sw  d  nw  sij,k − Vij,k s , = |qij,k | cCU,j − |qij,k | cij,k − Vij,k dt dt dt ij,k

(10)

12

where Vij,k , qij,k , and cij,k are the volume, volumetric flow rate, and concen-

13

nw −3 tration of k th edge of pore throat ij, respectively. ssw ij,k and sij,k [ML ] are

14

the concentrations of masses adsorbed at SW and NW interfaces, respec-

15

tively. 21

1

Within the pore throats we assume linear equilibrium adsorption, so that we

2

can write:

sw sw ssw ij,k = kd,ij,k aij,k cij,k nw nw snw ij,k = kd,ij,k aij,k cij,k ,

(11)

sw nw and kd,ij,k [L] are pore scale adsorption distribution coefficients where, kd,ij,k

at SW and NW interfaces, respectively. The specific surface area, aαw ij,k , is defined as, aαw ij,k

Aαw ij,k = , Vij,k

(12)

3

where Aαw ij,k is the total area of the appropriate αw interface (where α = s, n)

4

within the k th edge of pore throat ij.

5

Combining Eqs. (7) through (12) results in a linear set of equations to

6

be solved for cij,k and cCU,i . Since we discretize pore bodies and pore throats

7

on the basis of their saturation state (i.e., depending on whether a pore

8

is saturated or invaded by the nonwetting phase), the number of unknowns

9

depends on saturation. In the case of a fully saturated domain, the number of

10

unknowns is equal to Ntube + Nnode (Ntube is the total number of pore throats

11

and Nnode is the total number of pore bodies). In general, the number of pore

12

throats is larger than the number of pore bodies in pore network models. In

13

order to increase numerical efficiency, first, applying a fully implicit scheme,

14

we discretize Eq. (10) and determine cij,k in terms of cCU,i . This is then

15

substituted into the discretized form of Eq. (7), resulting in a set of equations

16

for cCU,i (Raoof and Hassanizadeh, 2010). This method considerably reduces

17

the computational times. For the accuracy of the scheme, the minimum time 22

1

step is chosen on the basis of residence times (Suchomel et al., 1998b), Δt ≤ min {Tij , Tij,k , TCU,i , Ti } ,

(13)

2

where T denotes the residence time pertaining to different elements within

3

the pore network. It is the volume of the element filled by the wetting phase

4

divided by the flux of wetting phase through the element.

5

After obtaining the solution at each time step, breakthrough curves (BTCs)

6

at a given longitudinal position are calculated by averaging the concentra-

7

tions of pore bodies centered at that position. In calculating BTCs, the

8

concentrations of pore bodies are weighted by their volumetric flow rate, re-

9

sulting in a flux-averaged concentration. That is, the normalized average

10

concentration, c(x, t), given by:  N x c(x, t) =

c (x, t)Q 1 i i i Ntx c0 i Qi t

i = 1, 2, 3, . . . , Nt ,

(14)

11

where c0 is inlet solute concentration and Ntx denotes the total number

12

of pore body elements that are located at the longitudinal coordinate x.

13

The longitudinal coordinate can be written as multiples of lattice size , i.e.

14

x = 1, 2, . . . , L. where  is the horizontal distance between centers of two

15

adjacent pore bodies. The breakthrough curve at the outlet is obtained by

16

plotting c(x = L, t).

17

3.3. Solute transport, including diffusion

18

For the case of reactive solute transport through the network including

19

both advection and diffusion, the mass balance Eq. (7) for a pore body

20

corner unit is extended to the following form: 23

ij

VCU,i dtd

(cCU,i ) =

throat N Nin edge  

j=1

ij

zi N edge  j=1 k=1

DAij,k

k=1

qij,k cij,k +

(cij,k −cCU,i ) lij

+

CU,i Nedge



n=1

CU,i Nin,edge



n=1

DAi,n

qi,n cCU,n − QCU,i cCU,i +

(cCU,n −cCU,i ) li,n

− RCU,i .

(15)

1

The fourth and fifth terms show the diffusive mass flux between corner

2

ij CU,i unit i and Nedge pore throats as well as Nedge corner units connected to it.

3

D is the molecular diffusion coefficient and Aij,k is the cross section of k th

4

edge of tube ij, which is a function of saturation. Ai,n is the cross section

5

of the pore body edge connecting corner unit i to corner unit n. We should

6

note that, for the case of saturated pores, there is no edge flow present and

7

the value of Nijedge will be equal to one.

8

Similarly, the mass balance Eq. (10) for an edge element of a drained

9

pore throat (assuming that corner unit j is the upstream node) is modified

10

to:

d (cij,k ) = |qij,k |cCU,j − |qij,k |cij,k + dt (cCU,j − cij,k ) (cCU,i − cij,k ) + DAij,k − Rij,k , DAij,k lij lij Vij,k

(16)

11

In Eqs. (15) and (16), RCU,i and Rij,k represent the reaction terms for corner

12

units of drained pore body i and k th edge of drained pore throat ij, respec-

13

tively. Note that for multi-component transport these parameters should be

14

l l rewritten as RCU,i and Rij,k for each chemical species l.

24

1

3.4. Biogeochemical Reaction Network Simulator (BRNS) PoreFlow uses Biogeochemical Reaction Network Simulator (BRNS) (Reg-

2 3

nier et al., 2002) for solving a set of coupled non-linear equations model-

4

ing reactive solutes. BRNS has been used as a flexible numerical engine

5

applicable for a variety of reactive transport problems in surface and sub-

6

surface environments such as, groundwater contamination, interaction be-

7

tween groundwater and sea water, seafloor carbon burial, and enhanced oil

8

recovery (Jourabchi et al., 2005; Spiteri et al., 2008b; Thullner et al., 2008;

9

Nick et al., 2012). One of the main features of BRNS is an automated proce-

10

dure for code generation using a MAPLE interface. This provides high flexi-

11

bility to include and combine alternative biogeochemical process descriptions

12

in the model. Using this interface, kinetically controlled and equilibrium re-

13

actions, and combinations of both can be specified with arbitrary form of

14

equations describing the transformations of the chemical species. A descrip-

15

tion on how symbolic programming can be used to create the model specific

16

source code necessary to the numerical solution of the governing equations

17

is explained in Regnier et al. (2002).

18

The reactive processes, rate parameters and equilibrium constants, species

19

involved, and stoichiometry are required to define a specific reaction network

20

and to solve for the reactive terms in each pore body and pore throats (i.e.,

21

reaction terms in equations 15 and 16). The reaction term c1 representing the

22

sum of reactions affecting concentration of a given species l can be written

23

as: c1

represents

25

l

R =

Nreaction 

al,m r l ,

(17)

m=1 1

where Nreaction is the number of reactions, r l represents, for kinetic reactions,

2

the rate of the reaction m and al,m is the stoichiometric coefficient of species l

3

in the reaction m. The form of rate expression, r l , is arbitrary, even nonlinear,

4

and can be a function of several concentrations within the system.

5

The MAPLE interface is employed to generate Fortran files, such as the

6

function residuals, which represent the reaction network, and the Jacobian

7

matrix, which contains the partial derivatives of the function residuals with

8

respect to the unknown concentrations. The generated Fortran files and the

9

solver engine for solving a nonlinear set of equations ensued from kinetic and

10

equilibrium reactions of multi components are embedded in PoreFlow. The

11

reactive solver uses a Newton-Raphson method to linearize the set of kinetic

12

and equilibrium reaction equations. In this work a sequential non-iterative

13

approach is implemented. Within the employed sequential approach, first

14

the transport step, and then the reactive step is solved in each time step of

15

the simulation.

16

4. Application examples

17

In this section, we present results for various simulations of (variably)

18

saturated flow and (reactive/adsorptive) transport using PoreFlow.

19

4.1. Non-reactive transport; saturated conditions

20

As the first example, we present results obtained using PoreFlow to sim-

21

ulate transport of a tracer solute through a column. We consider the case of 26

1

a pulse input of the solute. Averaging over the network cross section results

2

in a 1D concentration field. At the macro-scale, solute transport through a

3

column is modeled by the Advection Dispersion Equation (ADE), θ

∂2c ∂c ∂c + θv = θDL 2 , ∂t ∂x ∂x

(18)

4

where θ is porosity, v is average pore velocity, and DL is the longitudinal

5

dispersion coefficient. In this equation, porosity and v values are determined

6

directly from the pore network

7

the breakthrough curves calculated using PNM.

c1

and DL is estimated by fitting Eq. (18) c2 to

8

To verify the pore network model, we compare the resulting breakthrough

9

curve from the pore network model with the analytical solution of Eq. (18)

10

c3

representing the macro-scale behavior for the same value of porosity, aver-

11

age velocity, and the domain length as in the pore network. We use a network

12

of size Ni = 40, Nj = Nk = 20 (Model 1) with an average velocity of 0.72

13

m/d. The properties of the generated network are given in Table 1. A pulse

14

tracer concentration of one is injected through the inlet boundary for 1.8

15

pore volumes. At various time steps concentrations at the outlet face of the

16

network are averaged over the cross section to obtain the breakthrough curve.

17

Fig. 5 shows the resulting BTC together with the analytical solution. The

18

agreement with the analytical solution shows the accuracy of the numerical

19

scheme. c1 c2 c3

Text added. Text added. for the same value of porosity, average velocity, and the domain length as in the pore

network.

27

1

4.2. Non-reactive transport; variably saturated conditions An important parameter under variably saturated conditions is solute

2

c1

3

dispersion coefficient,

that depends not only on the Darcy velocity but also

4

on degree of saturation, which has been shown by Raoof and Hassanizadeh

5

(2013) c2 through pore scale modeling and comparison with experimental data

6

reported by Toride et al. (2003).

7

The degree of spreading of solutes is related to distribution of the pore

8

velocities within the pore spaces. Decrease in saturation creates corner flows

9

within drained pores and increases pore velocity fluctuations, which results

10

in a higher solute dispersion coefficient (Toride et al., 2003). Using Pore-

11

Flow, we can simulate transport of solutes under different saturation levels

12

(Raoof and Hassanizadeh, 2013). Since, resistance to the flow is included

13

within drained pore bodies and pore throats, different fluxes are calculated

14

for different edges of each drained pore, which are used in the calculation

15

of solute transport. Fig. 6 shows the BTCs of tracer solute obtained under

16

three different saturations for Model 1. There is a greater spreading of so-

17

lute (i.e. higher dispersivity) under lower saturations due to the presence of

18

longer flow paths, and wider distribution of arrival times. Such a behavior

19

has been observed through experiments (De Smedt et al., 1986; Toride et al.,

20

2003). c1

which depends not only on the Darcy velocity, but also on degree of saturation

c2

28

1

4.3. Adsorptive transport; variably saturated conditions

2

Adsorptive transport, such as transport of viruses and colloids, is of great

3

importance in studies of porous media (Bradford and Torkzaban, 2008). Un-

4

der variably saturated conditions, principal interactions usually occur not

5

only at the solid-wetting (SW) interfaces, but also at the nonwetting-wetting

6

(NW) interfaces (Torkzaban et al., 2008). These interactions are greatly

7

affected by the water content. Eq. (8) and (11) show that the concen-

8

tration of adsorbed solutes depends on specific interfacial areas, which, in

9

turn, are functions of saturation. Fig. 7 shows the normalized total inter-

10

facial areas as a function of saturation (Model 1), together with capillary

11

pressure-saturation as well as relative permeability-saturation curves. The

12

slope of kr −Sw is larger at higher saturation. This is due to the fact that the

13

nonwetting phase first occupies the larger pores (which have the most contri-

14

bution to the conductivity of the phase), after which there is a considerable

15

drop in conductivity of the wetting phase. Further decrease in saturation will

16

only cause slight decrease in permeability of the wetting phase, since almost

17

all larger and percolating pores are already drained. As shown in Fig. 7, the

18

SW interfacial area

19

the fluid films are neglected). This is due to the invasion of pores by the

20

nonwetting phase. Initially, under saturated conditions, there is no NW in-

21

terface present. Thus, the area of NW interface starts at zero, and increases

22

with decrease in saturation.

c3

decreases during drainage (interfaces associated with

At any given saturation, the fluids distribution within the pores and thus

23

c3

decrease

29

1

the amount of SW and NW interfaces in each pore are known. Utilizing this

2

information (e.g., using Eq. 12) the system of mass balance equations for

3

pore elements can be solved, and breakthrough curve can be obtained. Fig.

4

7 shows the BTC of an adsorbing solute,

5

the outlet of the pore network model with saturation of 0.5 .

6

4.4. Saturated reactive transport; seawater-fresh water mixing interface

c1

assuming kdsw = kdnw = 0.1mm, at

7

The reactive mixing between seawater and groundwater in coastal aquifers

8

influences the water quality of submarine groundwater discharge. While these

9

waters come into contact at the seawater groundwater interface by density

10

driven flow, their chemical components dilute and react through dispersion.

11

The interaction between intruding seawater and fresh groundwater causes

12

series of degradation processes in the coastal aquifer (Michael et al., 2003).

13

The multi-component reactive transport feature of PoreFlow permits to ad-

14

dress the chemical interactions that occur at seawater groundwater interface.

15

Here, we study these

16

intruding seawater and fresh groundwater at the pore scale level assuming

17

that the transverse dispersion is the main hydraulic process controlling such

18

mixing at this scale.

c1

interactions near a hypothetical interface between

19

To do so, a 3D pore network model of the size Ni = 301, Nj = Nk = 11 is

20

used (Model 2). The properties of the generated network are given in Table

21

1. The physical length of model domain is 300 mm, with a porosity value of

22

0.27 and permeability of 80 m/day. The model is assumed to be parallel to c1 c1

Text added. interaction

30

1

the freshwater-seawater interface along its length, with its lower boundary

2

representing the seawater. Following the work of Slomp and Van Cappellen

3

(2004) and Nick et al. (2012)c2 , we consider anoxic fresh water and oxic sea-

4

water conditions (Fig. 8)c3 . At the seawater boundary, seawater contains

5

O2 concentration of 0.2 mM and SO4 of 28 mM. The model is initially satu-

6

rated with fresh water. Then, water with the same density (1000 kg/m3 ) and

7

viscosity (1×10−3 Pa.s) but with different chemical compositions is injected

8

continuously through the left (inlet) boundary (Table 2). The dissolved or-

9

ganic carbon (DOC) of input water is (CH2 O)x · (NH3 )y · (H3 PO4 )z, where

10

x,y,z represent the C:N:P ratios, set to 106:11:1. Similar to Spiteri et al.

11

(2008a), we consider seven equilibrium/kinetic reactions which allows us to

12

study transformation processes affecting dissolved organic carbon degrada-

13

tion in saturated groundwater. These transformation processes are: aerobic

14

degradation, denitrification, nitrification, Fe(OH)3 reduction, SO2− reduc4

15

tion, Fe2+ oxidation and PO4 adsorption onto Fe(OH)3 (Table 3). Two dif-

16

ferent sets of simulations with average velocities of 2.317 and 0.464 m/yr

17

are performed. The top and bottom boundaries are set as no flow bound-

18

aries. Initially the concentrations of all the species are set to zero. The

19

aerobic degradation reaction rates is defined as given in Hunter et al. (1998)

20

and listed in Table 4. The limiting concentration for electron acceptors and

21

reaction rates are shown in Table 5.

22

units for immobile solutes are in mmol dm−3 , denoted as mM

23

Fig. 9 shows the average concentration profiles along the main flow path c2 c3

Reference added New figure added

31

1

for different components. The dissolved oxygen and sulphate diffuse laterally

2

(from the lower boundary representing the seawater) into the domain and

3

the reactive mixing zone width extends along the main flow path. Oxygen

4

concentration level is found to be lower for the simulation with lower flow

5

velocity. This is also reflected in breakthrough curves monitored at the outlet

6

(Fig. 10). As both examples are advective dominated flow, the advance of

7

tracer concentration through the lower boundary remains the same. This can

8

be seen in the tracer concentration profiles (Fig. 9 - O2 concentration profile)

9

and breakthrough curves for the tracer result (Fig. 10 - O2 breakthrough

10

curve). Therefore, under these conditions the molecular diffusion effect on

11

dispersion can be ignored.

12

As shown in Fig. 10, the oxygen concentration increases at early time

13

until the DOC plume propagates. Soon after the DOC plume reaches the

14

outlet boundary, the oxygen concentration decreases significantly. The sul-

15

fate plume travels nearly conservatively due to high concentration of this

16

component in the system. Denitrification is limited, because initially there is

17

no NO− 3 in the system and nitrification is relatively small and limited by the

18

oxygen supply. Fe2+ oxidation is also controlled by the available dissolved

19

oxygen. The DOC concentration reaching to the outlet for the case of lower

20

velocity is much less than that of the case with higher flow velocity. While

21

the velocity ratio between these two simulations is 5, the DOC mass flux

22

ratio is greater than 15. Therefore, the DOC degradation is more effective in

23

the model with lower velocity. Consequently N and P release rates are about

24

2 times larger for the higher velocity simulation. The NH4 produced in the

25

system due to primary redox reactions is nitrified to NO–3 when it comes in 32

1

contact with the O2 , which is dispersed laterally in the model. The removal

2

of NO–3 through denitrification is higher for the model with lower velocity.

3

The PO4 is removed through sorption on the Fe(OH)3 produced by the ox-

4

idation of Fe2+ where O2 is available. The removal of PO4 is enhanced for

5

the scenario with lower flow velocity because there is more mobile PO4 and

6

the Fe(OH)3 is available in the domain.

7

These results show that biochemical processes combined with flow and

8

transport can have a substantial effect on the release of nutrients into SGD.

9

Our main goal, however, is to illustrate the applicability of the PoreFlow for

10

modeling reactive transport in subsurface environment.

11

4.5. Reactive transport; including porosity change PoreFlow is capable of computing the change in porosity and the flow

12 13

c1

field due to the progress of chemical reactions (Raoof et al., 2012)c2 . It can

14

also be used to obtain the corresponding relationship between porosity and

15

permeability. In this section, we simulate a

16

including porosity changes. We consider dissolution of carbonate calcium,

17

CaCO3 , due to the injection of a low pH solution. This problem has various

18

applications such as CO2 sequestration into a reservoir rock.

c3

reaction/dissolution problem

19

The calcium carbonate equilibria parameters (Koutsoukos and Kontoyan-

20

nis, 1984) are given in Table 6. The dissolution of calcite has been proposed

21

to occur via three parallel reactions (Plummer et al., 1978; Chou et al., 1989) c1 c2 c3

filed Reference added recation

33

2+ − −− CaCO3 (s) + H+  − − Ca + HCO3

(19)

2+ − −− CaCO3 (s) + H2 CO3  − − Ca + 2 HCO3

(20)

2+ 2− −− CaCo3 (s) +  − − Ca + CO3

(21)

1

The reaction consumes H+ , and releases Ca2+ and carbonate. The rate

2

of calcite dissolution is described by a Transition State Theory (TST) rate

3

law (Chou et al., 1989): Rate = A(k1 aH+ + k2 aH2 CO3

  IAP + k3 ) 1 − , Keq

(22)

4

where k1 , k2 , and k3 are the reaction rate constants (Li et al., 2008), aH+ and

5

aH2 CO3 are the activities of hydrogen ion and carbonic acid, IAP is the ion

6

activity product defined as aCa2+ aCO2− , Keq is the equilibrium constant for 3

7

the reaction (Chou et al., 1989), and A is the reactive surface area. Eq. (22)

8

is used to calculate the dissolution rate within a pore.

9

We use concentrations in equilibrium with CaCO3 as the initial condi-

10

tions. The diffusion coefficient is set to 1.0×10−9 [m2 /s] and average velocity

11

is 0.1 [m/d]. Due to dissolution of calcite, the initial (equilibrium) pH is as

12

high as 10. Subsequent injection of solute with pH of 3.0 into a network of

13

size Ni = Nj = Nk = 10 (Model 3) for 60 days caused the dissolution of

14

calcite, and therefore, an increase in porosity.

15

radius of the pore throats (that affect permeability) and pore bodies (which

16

affect the storage of chemical products) is changed at each time step due

17

to the dissolution-precipitation processes. Having the dissolved/precipitated c1

Text added.

34

c1

During the simulation, the

1

mass of solid-phases calculated at each time step, the change in pore vol-

2

umes is calculated for each pore using density of each solid phase.Then, the

3

radius of each pore is updated accordingly. The above procedure is repeated

4

for each time increment to calculate the evolution of cement permeability

5

during the degradation process. By averaging over the network domain at

6

successive time steps, we calculate the evolution of porosity and permeabil-

7

ity as chemical reactions progress. Dissolution of calcite causes an increase

8

in pore throat sizes which have the most effect on permeability of the pore

9

network. The initial pore throat size distribution together with the pore

10

throat size distribution after dissolution simulation are shown in Fig. 11a.

11

Fig. 11b shows the corresponding increase in permeability due to evolution

12

in porosity. Results show a five-fold increase in permeability as a result of

13

5.0% increase in porosity.

14

5. Conclusion

15

In this study, we introduce a pore-scale modeling tool, PoreFlow, to quan-

16

tify and describe the physical/chemical processes that govern the transport

17

of both tracer and reactive/adsorptive solutes in (variably-)saturated porous

18

media. PoreFlow is a FORTRAN-90 modular package whose capabilities in-

19

clude the generation of random structure or user defined networks, simulation

20

of drainage process, discretization of pore spaces on the basis of saturation

21

state of each pore, and solving flow and transport of multi-component reac-

22

tive solutes under saturated and variably saturated conditions, using several

23

complex algorithms. The solute transport governing equations are solved

24

applying a fully implicit numerical scheme. An efficient substitution meth35

1

ods is applied to make the algorithm computationally more effective and

2

appropriate for parallel computations. In presence of the nonwetting phase,

3

there are fluid-fluid (i.e., nonwetting-wetting) interfaces in addition to solid-

4

fluid (i.e., solid-wetting) interfaces. So, adsorption to fluid-fluid interfaces

5

may occur in addition to adsorption to solid-fluid interfaces. Coupling with

6

Biogeochemical Reaction Network Simulator (BRNS), make it possible to

7

simulate complex system of equilibrium or kinetic reactions. The reactions

8

can be heterogenous resulting in evolution of porosity and permeability of

9

the porous structure.

10

We present results of various simulations showing the capability of Pore-

11

Flow to model reactive flow and transport in porous media. This helps un-

12

derstanding of the processes, to perform upscaling, and to develop macroscale

13

constitutive relations suitable as inputs for continuum scale models. These

14

include the relationships between phase saturation and parameters such as

15

capillary pressure (Pc − Sw curve), interfacial area (a − Sw curve), relative

16

permeability (kr − Sw curve) and solute dispersivity (α − Sw curve), as well

17

as the relationship between porosity and permeability (θ − k curve), and

18

tortuosity or effective diffusion and porosity (Def f − θ curve).

36

1 2 3

Acharya, R., van der Zee, S., and Leijnse, A. (2005). Transport modeling of nonlinearly adsorbing solutes in physically heterogeneous pore networks. Water Resources Research, 41(2):W02020.

7

Algive, L., B´ekri, S., Nader, F. H., Lerat, O., and Vizika, O. (2012). Impact of diagenetic alterations on the petrophysical and multiphase flow properties of carbonate rocks using a reactive pore network modeling approach. Oil & Gas Science and Technology–Revue dIFP Energies nouvelles, 67(1):147–160.

8

Bear, J. (1972). Dynamics of fluids in porous media, 764 pp.

9

Bear, J. (1988). Dynamics of fluids in porous media. Dover Publications.

4 5 6

10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34

Bijeljic, B., Mostaghimi, P., and Blunt, M. (2011). Signature of non-fickian solute transport in complex heterogeneous porous media. Physical Review Letters, 107(20):204502. Blunt, M. (2001). Flow in porous media; pore-network models and multiphase flow. Current Opinion in Colloid & Interface Science, 6(3):197–207. Bradford, S. and Torkzaban, S. (2008). Colloid transport and retention in unsaturated porous media: A review of interface-, collector-, and pore-scale processes and models. Vadose Zone Journal, 7(2):667–681. Brooks, R. H. and Corey, A. T. (1964). Hydraulic properties of porous media. Hydrol. Pap. 3, Colorado State University. Celia, M., Reeves, P., and Ferrand, L. (1995). Recent advances in pore scale models for multiphase flow in porous media. Reviews of Geophysics, 33(2):1049. Chandler, R., Lerman, K., Koplik, J., and Willemsen, J. (1982). Capillary displacement and percolation in porous media. Journal of Fluid Mechanics, 119:249–67. Chatzis, I. and Dullien, F. (1985). The modeling of mercury porosimetry and the relative permeability of mercury in sandstones using percolation theory. Int. Chem. Eng.;(United States), 25(1). Chou, L., Garrels, R., and Wollast, R. (1989). Comparative study of the kinetics and mechanisms of dissolution of carbonate minerals. Chemical Geology, 78(3-4):269–282. Class, H., Helmig, R., and Neuweiler, I. (2008). Sequential coupling of models for contaminant spreading in the vadose zone. Vadose Zone Journal, 7(2):721. De Smedt, F., Wauters, F., and Sevilla, J. (1986). Study of tracer movement through unsaturated sand. Geoderma, 38(1-4):223–236. Dupin, H., Kitanidis, P., and McCarty, P. (2001). Pore-scale modeling of biological clogging due to aggregate expansion: a material mechanics approach. Water Resources Research, 37(12):2965–2979.

37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35

Fatt, I. (1956). The network model of porous media. Trans. Am. Inst. Mem. Metall Pet. Eng, 207:144–181. Geiger, S., Schmid, K., and Zaretskiy, Y. (2012). Mathematical analysis and numerical simulation of multi-phase multi-component flow in heterogeneous porous media. Current Opinion in Colloid & Interface Science. Gharasoo, M., Centler, F., Regnier, P., Harms, H., and Thullner, M. (2011). A reactive transport modeling approach to simulate biogeochemical processes in pore structures with pore-scale heterogeneities. Environmental Modelling & Software. Held, R. and Celia, M. (2001). Pore-scale modeling and upscaling of nonaqueous phase liquid mass transfer. Water Resources Research, 37(3):539–549. Hunter, K., Wang, Y., and Van Cappellen, P. (1998). Kinetic modeling of microbiallydriven redox chemistry of subsurface environments: coupling transport, microbial metabolism and geochemistry. Journal of Hydrology, 209(1-4):53 – 80. Ioannidis, M., Kwiecien, M., Chatzis, I., MacDonald, I., and Dullien, F. (1997). Comprehensive pore structure characterization using 3d computer reconstruction and stochastic modeling. In Society of Petroleum Engineers; Annual Technical Conference and Exhibition, page 9. Jivkov, A. P., Hollis, C., Etiese, F., McDonald, S. A., and Withers, P. J. (2013). A novel architecture for pore network modelling with applications to permeability of porous media. Journal of Hydrology. Jourabchi, P., Regnier, P., and Van Cappellen, P. (2005). Quantitative interpretation of ph distributions in aquatic sediments; a reaction-transport modeling approach. American Journal of Science, 305(9):919–956. Kang, Q., Lichtner, P. C., Viswanathan, H. S., and Abdel-Fattah, A. I. (2010). Pore scale modeling of reactive transport involved in geologic co2 sequestration. Transport in porous media, 82(1):197–213. Kechagia, P., Tsimpanogiannis, I., Yortsos, Y., and Lichtner, P. (2002). On the upscaling of reaction-transport processes in porous media with fast or finite kinetics. Chemical engineering science, 57(13):2565–2577. Kim, D. and Fogler, H. (2000). Biomass evolution in porous media and its effects on permeability under starvation conditions. Biotechnology and bioengineering, 69(1):47– 56. Kim, D., Peters, C., and Lindquist, W. (2011). Upscaling geochemical reaction rates accompanying acidic co2-saturated brine flow in sandstone aquifers. Water Resources Research, 47(1).

38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35

Knutson, C., Werth, C., and Valocchi, A. (2001). Pore-scale modeling of dissolution from variably distributed nonaqueous phase liquid blobs. Water Resources Research, 37(12). K¨ ohne, J., Schl¨ uter, S., and Vogel, H. (2011). Predicting solute transport in structured soil using pore network models. Vadose Zone Journal, 10(3):1082–1096. Koutsoukos, P. and Kontoyannis, C. (1984). Precipitation of calcium carbonate in aqueous solutions. J. Chem. Soc., Faraday Trans. 1, 80(5):1181–1192. Li, L., Peters, C., and Celia, M. (2006). Upscaling geochemical reaction rates using porescale network modeling. Advances in Water Resources, 29(9):1351–1370. Li, L., Steefel, C., and Yang, L. (2008). Scale dependence of mineral dissolution rates within single pores and fractures. Geochimica et Cosmochimica Acta, 72(2):360–377. Lichtner, P. (1985). Continuum model for simultaneous chemical reactions and mass transport in hydrothermal systems. Geochimica et Cosmochimica Acta. Lindquist, W., Venkatarangan, A., Dunsmuir, J., and Wong, T. (2000). Pore and throat size distributions measured from synchrotron x-ray tomographic images of fontainebleau sandstones. JOURNAL OF GEOPHYSICAL RESEARCH-ALL SERIES-, 105(B9):21– 21. Mason, G. and Morrow, N. (1991). Capillary behavior of a perfectly wetting liquid in irregular triangular tubes. J. Colloid Interface Sci, 141(1):262–274. Mehmani, Y., Sun, T., Balhoff, M., Eichhubl, P., and Bryant, S. (2012). Multiblock porescale modeling and upscaling of reactive transport: Application to carbon sequestration. Transport in Porous Media, 95(2):305–326. Michael, H., Lubetsky, J., and Harvey, C. (2003). Characterizing submarine groundwater discharge: A seepage meter study in waquoit bay, massachusetts. Geophysical Research Letters, 30(6). Nick, H. and Matth¨ ai, S. (2011). Comparison of three FE-FV numerical schemes for singleand two-phase flow simulation of fractured porous media. Transport in Porous Media, pages 1–24. Nick, H., Paluszny, A., Blunt, M., and Matthai, S. (2011). Role of geomechanically grown fractures on dispersive transport in heterogeneous geological formations. Physical Review E, 84(5):056301. Nick, H., Raoof, A., Centler, F., Thullner, M., and Regnier, P. (2012). Reactive dispersive contaminant transport in coastal aquifers: Numerical simulation of a reactive henry problem. Journal of contaminant hydrology. Oren, P. and Bakke, S. (2003). Reconstruction of berea sandstone and pore-scale modelling of wettability effects. Journal of Petroleum Science and Engineering, 39(3-4):177–199.

39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35

Oren, P., Bakke, S., and Arntzen, O. (1998). Extending predictive capabilities to network models. SPE Journal, 3(4):324–336. Ovaysi, S. and Piri, M. (2011). Pore-scale modeling of dispersion in disordered porous media. Journal of contaminant hydrology. Plummer, L., Wigley, T., and Parkhurst, D. (1978). The kinetics of calcite dissolution in c 02-water systems at 5 deg to 60 deg c and 0. 0 to 1. 0 atm c 02. American Journal of Science, 278. Raoof, A. (2011). Reactive/adsorptive transport in (partially-) saturated porous media; from pore scale to core scale. Utrecht University. Raoof, A. and Hassanizadeh, S. (2009). A new method for generating pore-network models of porous media. Transport in Porous Media, Volume 81, Number 3:391–407. Raoof, A. and Hassanizadeh, S. (2010). Upscaling transport of adsorbing solutes in porous media: Pore-network modeling. Vadose Zone Journal, Accepted for publication. Raoof, A. and Hassanizadeh, S. (2011). A new formulation for pore-network modeling of two-phase flow. WRR. Raoof, A. and Hassanizadeh, S. (2013). Saturation dependent solute dispersivity in porous media; pore-scale processes. Water Resources Research. Raoof, A., Nick, H., Wolterbeek, T., and Spiers, C. (2012). Pore-scale modeling of reactive transport in wellbore cement under co2 storage conditions. International Journal of Greenhouse Gas Control. Regnier, P., O’Kane, J., Steefel, C., and Vanderborght, J. (2002). Modeling complex multi-component reactive-transport systems: towards a simulation environment based on the concept of a knowledge base. Applied Mathematical Modelling, 26(9):913 – 927. Rosenzweig, R., Furman, A., and Shavit, U. (2013). A channel network model as a framework for characterizing variably saturated flow in biofilm-affected soils. Vadose Zone Journal, 12(2). Sholokhova, Y., Kim, D., and Brent Lindquist, W. (2009). Network flow modeling via lattice-boltzmann based channel conductance. Advances in Water Resources, 32(2):205– 212. Slomp, C. and Van Cappellen, P. (2004). Nutrient inputs to the coastal ocean through submarine groundwater discharge: controls and potential impact. Journal of Hydrology, 295(1-4):64–86. Spiteri, C., Slomp, C., Charette, M., Tuncay, K., and Meile, C. (2008a). Flow and nutrient dynamics in a subterranean estuary (waquoit bay, ma, usa): Field data and reactive transport modeling. Geochimica et Cosmochimica Acta, 72(14):3398 – 3412.

40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31

Spiteri, C., Slomp, C., Tuncay, K., and Meile, C. (2008b). Modeling biogeochemical processes in subterranean estuaries: the effect of flow dynamics and redox conditions on submarine groundwater discharge. Water Resources Research, 44(2):18 pp. Suchomel, B., Chen, B., and Allen, M. (1998a). Network model of flow, transport and biofilm effects in porous media. Transport in porous media, 30(1):1–23. Suchomel, B., Chen, B., and Allen, M. (1998b). Network model of flow, transport and biofilm effects in porous media. Transport in porous media, 30(1):1–23. Thullner, M., Kampara, M., Richnow, H. H., Harms, H., and Wick, L. Y. (2008). Impact of bioavailability restrictions on microbially induced stable isotope fractionation. 1. theoretical calculation. Environmental Science & Technology, 42(17):6544–6551. Toride, N., Inoue, M., and Leij, F. (2003). Hydrodynamic dispersion in an unsaturated dune sand. Soil Science Society of America Journal, 67(3):703–712. Torkzaban, S., Bradford, S., van Genuchten, M., and Walker, S. (2008). Colloid transport in unsaturated porous media: The role of water content and ionic strength on particle straining. Journal of contaminant hydrology, 96(1-4):113–127. Van Genuchten, M. (1980). A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Science Society of America Journal, 44(5):892–898. Varloteaux, C., B´ekri, S., and Adler, P. M. (2012). Pore network modelling to determine the transport properties in presence of a reactive fluid: from pore to reservoir scale. Advances in Water Resources. Vasilyev, L., Raoof, A., and Nordbotten, J. M. (2012). Effect of mean network coordination number on dispersivity characteristics. Transport in Porous Media, 95(2):447–463. Vogel, H. and Roth, K. (2001). Quantitative morphology and network representation of soil pore structure. Advances in Water Resources, 24(3-4):233–242. Wilkinson, D. and Willemsen, J. (1983). Invasion percolation: a new form of percolation theory. Journal of Physics A: Mathematical and General, 16:3365–3376. Yu, L. and Wardlaw, C. (1986). Mechanisms of nonwetting phase trapping during imbibition at slow rates. Journal of Colloid and Interface Science, 109(2):473–486. Zhou, D., Dillard, L., and Blunt, M. (2000). A physically based model of dissolution of nonaqueous phase liquids in the saturated zone. Transport in Porous Media, 39(2):227– 255.

41

42

Figure 1: Program flowchart and structure in PoreFlow.

Figure 2: Schematic of a network consisting of three pore bodies in each direction (i.e., Ni=3, Nj=3, Nk=3). Numbers inside squares denote all possible throat directions and plain numbers are pore body numbers. To keep the figure less crowded, only throats which are connecting pore body 14 to its neighboring pore bodies are shown [After (Raoof and Hassanizadeh, 2009)].

Probability

0.2

a

0.15 0.1 0.05 0

1

3

5 7 9 11 13 Coordination number

15

Figure 3: The coordination number distribution and a representative domain of MultiDirectional Pore-Network (MDPN).

c2

The color-bar shows the coordination number.

Table 1: Properties of three pore networks produced by MDPN.

43

PoreFlow-Tables.tex

Characteristic Initial Nthroat (before elimination) Initial Nnode (before elimination) Nthroat (after elimination) Nnode (after elimination) Connections to the inlet boundary Connections to the outlet boundary pore bodies on the inlet boundary pore bodies on the outlet boundary Porosity

Table 1: Model 1 187512

Model 2 336332

Model 3 9792

16000

30000

1000

33239 14804 578

86888 29257 215

2931 977 237

581

193

265

315

91

93

313

86

92

0.25

0.27

0.12

Table 2: DOC 0.0 1.0

Seawater groundwater interface boundary Inland groundwater plume

1

NO3– 0.0 0.0

NH+ 4 0.0 0.0

PO4 0.0 0.0

O2 0.05 0.0

Fe 2+ 0.0 0.1

SO42 – 28 0.0

Figure 4: (a) Example of two drained pore-bodies connected to each other using a drained pore throat with angular cross section. Each drained pore body is refined into 8 corner units (pore elements), and (b) a pore body which is invaded by the nonwetting phase (dark grey) through one of its throats, and, as a result, reduces the connectivity of the neighboring saturated pore throats. The light grey shade represents the wetting fluid phase.

1 Network model Analytical solution

Concentration [−]

0.8

0.6

0.4

0.2

0 0

1

2

3

4

Pore volume

Figure 5: Comparison of BTC obtained from pore network (Model 1) with analytical solution for a pulse input of solute to a column. The solid line shows average concentration at the outlet of the network and the symbols are analytical solution of continuum scale Eq. (18).

44

1

S = 1.0 w

S = 0.8 w

0.8

S = 0.5

Concentration

w

0.6

0.4

0.2

0 0

1

2

3

4

5

6

Pore volume

Figure 6: BTCs obtained from the pore network (Model 1) at different saturations: Sw = 1.0, 0.8, and 0.55. BTCs obtained under lower saturations show more spreading compared with BTC of tracer for saturated condition.

Table 2: Boundary concentrations at the left and bottom boundaries of Model 2, units for mobile solutes are in mmol dm−3 pore water, denoted as mM. DOC

NO–3

NH+ 4

PO4

O2

Fe2+

SO2– 4

Seawater groundwater interface boundary

0.0

0.0

0.0

0.0

0.05

0.0

28

Inland groundwater plume

1.0

0.0

0.0

0.0

0.0

0.1

0.0

45

700

1

a

500

400

300

200

b

0.8

Relative permeability

Capillary pressure [Pa]

600

0.6

0.4

0.2

0

0.2

0.4

0.6

0.8

0 0

1

0.2

0.4

Saturation 1

1

SW interface NW interface

0.6

0.4

0.2

0 0

0.8

1

d

0.8

c

Concentration [−]

Normalized interface

0.8

0.6

Saturation

0.6

0.4

0.2

0.2

0.4

0.6

0.8

0 0

1

10

Saturation

20

30

40

50

60

70

Pore volume

Figure 7: (a) Capillary pressure-saturation curve, (b) relative permeability-saturation curve, (c) the normalized total interfacial area of solid-wetting phases, SW, and nonwetting-wetting phases, NW, as a function of average wetting saturation, Sw (the areas are normalized to the total area under saturated conditions, which only belongs to the solid-wetting surfaces), and (d) the resulting BTC from the network (Model 1) under variably saturated conditions. The BTC belong to kdsw = kdnw = 0.1mm and Sw = 0.5. Interfaces associated with the fluid films are neglected.

46

Tracer concentration

boundary No flow

1

0.9 0.8

Anoxic groundwater discharge

0.7 0.6 0.5 0.4 0.3 0.2 0.1

undary -

c2

ace bo ter interf

w

oxic sea

wa

r ground

Seawate

Figure 8:

)

oundary

flow b ater (No

A sample of the pore network domain employed for seawater-fresh water

problem. The concentration front represents a continuous tracer injection from the left boundary. The concentrations at the model boundaries are listed in Table 2.

47

0

Table 3: The reaction network for Model 2. Primary redox reactions Aerobic degradation (CH2 O)106 ·(NH3 )11 ·(H3 PO4 ) + 106 O2 −−→ 2– 9 HCO–3 + 97CO2 + 11 NH+ 4 + HPO4 + 97 H2 O

Denitrification (CH2 O)106 ·(NH3 )11 ·(H3 PO4 ) + 84.8 NO–3 −−→ 2– 42.4 N2 + 12.2 CO2 + 93.8 HCO–3 + 11 NH+ 4 + HPO4 + 54.6 H2 O

Fe(OH)3 reduction (CH2 O)106 ·(NH3 )11 ·(H3 PO4 ) + 424 Fe(OH)3 + 751 CO2 −−→ 2– 424 Fe2+ + 857 HCO–3 + 11 NH+ 4 + HPO4 + 309 H2 O

Sulfate reduction (CH2 O)106 ·(NH3 )11 ·(H3 PO4 ) + 53 SO2– 4 + 9 CO2 + 9 H2 O −−→ 2– 53 H2 S + 115 HCO–3 + 11 NH+ 4 + HPO4

Secondary redox reactions Nitrification – – NH+ 4 + 2 O2 + 2 HCO3 −−→ NO3 + 2 CO2 + 3 H2 O

Fe2+ oxidation Fe2+ + 2 HCO–3 + 0.25 O2 + 0.5 H2 O −−→ Fe(OH)3 + 2 CO2 Equilibrium sorption P adsorption PO4aq ←−→ PO4ads

48

Table 4: The kinetic/equilibrium formulations for Model 2. Reaction

Condition

Rates

Aerobic degradation

If O2 > KmO2

KDOC [DOC]

else

KDOC [DOC]

If O2 > KmO2

0

else If NO–3 > KmNO−

KDOC [DOC]

else If NO–3 < KmNO−

KDOC [DOC]

If NO–3 > KmNO−

0

else If Fe(OH)3 > KmFe(OH)

KDOC [DOC]

Denitrification

3

3

Fe(OH)3 reduction

3

3

Sulfate reduction

If Fe(OH)3 > KmFe(OH)

0

3

else If else If

SO2– 4

> KmSO2− < KmSO2−

1−

[O2 ] [KmO ]





KDOC [DOC]

4

[O2 ] [KmO ] 2

 

2

[NO− 3 ] [Km −] NO3

1−

[O2 ] [KmO ]

1−

[O2 ] [KmO ]

2

2

−

[NO− 3 ] [Km −]

−

[NO− 3 ] [Km −]



KDOC [DOC]

4



1−



KDOC [DOC]

SO2– 4

2



else If Fe(OH)3 < KmFe(OH)

3

[O2 ] [KmO ]

1− 1−

Nitrification

Knit [NH4+][O2]

Fe2+ oxidation

KFe [Fe2+ ][O2]

P adsorption

Kd =

[O2 ] [KmO ] 2

[O2 ] [KmO ] 2

− −

NO

NO

3

  [Fe(OH)3 ] [KmFe(OH) ] 3

3

[NO− 3 ] [Km −] NO3

[NO− 3 ] [Km −] NO 3

 [Fe(OH)3 ] [KmFe(OH) ]

−

3

−

[Fe(OH)3 ] [KmFe(OH) ] 3

[Pads ] [Paq ]

Table 5: Parameter values used for simulation of the DOC degradation. KDOC is the rate constant for decomposition of (CH2 O)106 ·(NH3 )11 ·(H3 PO4 ); Knit is the rate constant of nitrification; KFe is the rate for Fe2+ reoxidation; KmO2 , KmNO− , KmFe(OH)3 and 3

KmSO2− are limiting concentrations of O2 , NO–3 , Fe(OH)3 and SO2– 4 respectively; Kd is 4

the adsorption coefficient (Hunter et al., 1998). KDOC

Knit

KFe

s−1

mM−1 s−1

s−1

KmO2

KmNO−

mM

mM

mmol

1.0×10−7

1.0×10−4

1.0×10−4

3.0×10−2

1.0×10−3

20

3

49

KmFe(OH)3 dm−3

KmSO2−

Kd

mM

dm3 mmol−1

1.6

100

4



[SO2− 4 ] [Km 2− ] SO

4

b

a

c

c

Figure 9: a) DOC, b) O2 , c) PO4 aq , and d)NO− 3 concentration profiles along the main flow path plotted after 1.86 pore volume. Units for mobile solutes are in mmol dm−3 pore volume. Dashed lines are the profiles for the model with the velocity of 2.317 m/yr and the solid lines represents the result of the model with velocity value of 0.464 m/yr. The lines marked with the circular symbols represent the result for the non reactive tracer.

50

a

b

c

d

2+ Figure 10: Flux weighted breakthrough curves for a) DOC, b) O2 , c) SO2− are 4 , and d)Fe

illustrated for different flow velocities after 1.86 pore volume. The dashed lines indicate the result for the model with the velocity of 2.317 m/yr and the solid line represents the result of model with velocity value of 0.464 m/yr. The lines marked with the circular symbols represent the result for the non-reactive tracer. Units for mobile solute concentrations (vertical axes) are in mmol dm−3 pore water (mM).

51

Table 6: Calcium carbonate equilibria at 25◦ C (Koutsoukos and Kontoyannis, 1984) K

equilibrium −

H2 O − − − − H+ + OH

1.007×10−14

3 −− H+ + HCO− 3 −− H2 CO

2.249 × 106

− H+ + CO2− − − − HCO3 3 −

2.133 × 1010

− −− + Ca+ 2 + HCO3 −− CaHCO3

10.0

Ca+ 2

−

+ OH − − − − CaOH

+

19.95

−− CO2 (aq)  − − H2 CO3

3.390 × 10−2

1

0.3

5.5

a

5 4.5

Conductivity change

Normalized frequency

0.25

b

0.2 0.15 0.1

4 3.5 3 2.5 2

0.05 1.5 0 0

0.02

0.04

0.06

0.08

1 0.12

0.1

0.13

Pore throat radius [mm]

0.14

0.15

0.16

0.17

0.18

Porosity

Figure 11: a) Initial distribution of pore throats together with pore throat size distribution after the dissolution reaction, b) relation between permeability increase with porosity. k0 = 1.2 × 10−12 m2 is the initial permeability for Model 3.

52