From pore scale to continuum scale modeling of infiltration

Accepted Manuscript

From pore scale to continuum scale modeling of infiltration J. Tzavaras, M. Kohne, H.-J. Vogel ¨ PII: DOI: Reference:

S0309-1708(16)30281-0 10.1016/j.advwatres.2017.03.005 ADWR 2795

To appear in:

Advances in Water Resources

Received date: Revised date: Accepted date:

23 July 2016 15 November 2016 7 March 2017

Please cite this article as: J. Tzavaras, M. Kohne, H.-J. Vogel, From pore scale to continuum scale ¨ modeling of infiltration, Advances in Water Resources (2017), doi: 10.1016/j.advwatres.2017.03.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 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 • pore network model based on cylindrical pores and nodes without spatial extent • implementation of an infiltration process

AC

CE

PT

ED

M

AN US

• comparison to Richards equation (continuum scale modeling)

CR IP T

• biconjugate gradient method with pressure field iterations

1

ACCEPTED MANUSCRIPT

From pore scale to continuum scale modeling of infiltration J. Tzavarasa,∗, M. K¨ohnea , H.-J. Vogela Centre for Environmental Research, Theodor-Lieser-Str. 4, 06120 Halle (Saale), Germany

CR IP T

a UFZ-Helmholtz

Abstract

AN US

Water infiltration in soil at the continuum scale is typically modeled using Richards equation. This requires effective material properties, the water retention characteristic and the unsaturated hydraulic conductivity function. During infiltration the gaseous phase is replaced by water within a complex porous structure. This produces phenomena such as irregular infiltration fronts, hydraulic non-equilibrium and hysteresis which are ultimately related to pore scale processes. In this work we simulate infiltration at the pore scale using a pore network model where the pore structure can be adapted to real structures in terms of pore size distribution and pore topology. We compare the horizontally averaged dynamics of water content and water potential to the results obtained from Richards’ equation. This provides evidence that the pore network representation is consistent with the typical macroscopic dynamics at the continuum scale. However, assuming immiscible and incompressible fluid phases leads to unrealistic gas entrapment in the network model. We conclude that mechanisms leading to a release of trapped gas phase in natural porous media are difficult to represent in a pore network model. Yet, the proposed pore scale approach has the potential to study fundamental relations between pore structure and macroscopic phenomena not only for drainage but also for infiltration.

2

1. Introduction

ED

1

M

Keywords: pore network model, pore scale modeling, continuum scale modeling, infiltration, Richards equation

Pore network models (PNMs) have long been used to describe water dynamics in porous media. Based

5

Joekar-Niasar and Hassanizadeh [4]. Two-phase flow in a steady state regime was simulated by Constan-

6

tinides and Payatakes [5] to predict the macroscopic relative permeabilities, and to relate relative permeabil-

PT

4

on the pioneering work of Fatt [1, 2, 3] there is an enormous body of literature mostly dealing with drainage processes (i.e. non-wetting phase replacing wetting phase) and stationary flow. A recent review is given by

3

ities to the flow mechanisms at pore level. Wetting process are much more difficult to describe, since the

8

dynamics of any individual interface need to be modeled. This can be done by pore scale simulations based

9

on Stokes or Navier-Stokes equation to be solved numerically based on Lattice Boltzmann (Ahrenholz et al.

10

[6], Porter et al. [7]) or similar approaches. The major drawback is that the pore space need to be explic-

AC

CE

7

11

itly represented within a very limited volume and at a limited resolution. This is why the approach based ∗ Corresponding

author Email addresses: [email protected] (J. Tzavaras), [email protected] (M. K¨ohne), [email protected] (H.-J. Vogel)

Preprint submitted to Elsevier

July 2016

ACCEPTED MANUSCRIPT

12

on Navier-Stokes equation is typically applied to mono-scale porous media as sand or sandstone where the

13

relevant pore structure can be captured within a few cubic millimeters. To overcome this problem, Yang

14

et al. [8] added a Darcy-term to represent flow in the porous matrix. In this way, smaller pores that are not

17

wetting process is difficult to compute. Wetting processes, especially infiltration into relatively dry soil leads

18

to macroscopic phenomena as hydraulic non-equilibrium or hysteresis which are relevant for the redistribu-

19

tion of precipitation water in soil and which are deemed to be explained by effects of heterogeneities at the

20

pore scale. In petroleum engineering this is a highly relevant problem for replacing oil by water.

CR IP T

16

explicitly represented for solving Navier-Stokes equation could be represented. A pore network model as used in our study offers the possibility to represent more complex structures as typical for soils, however the

15

In contrast to static drainage processes, the simulation of infiltration into network models is challenging

22

since the pressure distribution has to be recalculated after each change in fluid saturation to correctly repre-

23

sent the final fluid configuration at equilibrium. Solution to this problem has been provided by Mogensen

24

and Stenby [9] and was more thoroughly described and elaborated in Mogensen et al. [10] who introduced

25

an iterative method to recalculate the pressure field after each time step of infiltration. More recently another

26

approach was presented by Tørå et al. [11] who expanded the solution for two-phase flow in a pore network

27

model to imbibition, drainage and steady state displacement.

AN US

21

In our paper we apply the concept of Mogensen et al. [10] to water infiltration in soil. Our pore network

29

model (PNM) is based on straight cylindrical pores connected to nodes without spatial extent. The main

30

motivation is to explore if the pore scale description is consistent with continuum scale models in terms of

31

the macroscopic dynamics of hydraulic state variables, i.e. water content and water potential.

32

2. Pore scale modeling

33

2.1. Network generation

ED

M

28

To generate a realistic pore network model (PNM) we use the concept introduced by Vogel and Roth

35

[12]. Thereby, relevant structural properties as pore size distribution and pore topology can be preserved

36

while less relevant details such as surface roughness are simplified. The PNM is based on a face-centered

37

cubic grid (coordination number 12) with cylindrical bonds representing pores while the nodes have zero

38

volume. Pore size distribution and pore topology are adjusted to measured values as obtained from 3D

39

CT images. To represent the measured pore geometry only a small part of the 12 possible bonds per node

40

are required so that the effective coordination number Zeff is smaller than Zmax = 12. In contrast to other

41

approaches where bonds are removed by some random or regular process (Raoof and Hassanizadeh [13])

42

we build up the network starting from the largest pores so that the Euler number of each pore size class as

43

measured in the CT images is reproduced. A special feature of our network model is that the remaining

44

bonds (i.e. Zmax − Zeff ) are considered as “background porosity”. This accounts for the fact that in soil

AC

CE

PT

34

45

there are always small pores that are not resolved by the imaging technique (CT). The resolution of CT

46

images is at best 0.01 mm while clay particles in soil are much smaller. As a consequence, these small pores

47

cannot be represented in the network of discrete pores, however, they definitely contribute to water flow. 2

ACCEPTED MANUSCRIPT

pore diameter EPC (str.) EPC (ran.)

0.04 -34.49 -45.37

0.066 -1.92 -5.61

0.092 -0.01 1.57

0.118 0.47 3.69

0.144 0.37 2.65

0.169 0.39 1.66

0.207 0.27 1.05

0.253 0.24 0.88

0.299 0.15 0.67

0.345 0.08 0.54

Table 1: Pore diameter and measured Euler Numbers for the structured and random PNM.

48

With adding this background porosity described by some mean pore size smaller than the other pores in the PNM, this effect is accounted for. In this way also the water phase is always continuous reflecting the natural

50

situation of moist soil. In other PNMs this important feature is obtained by using angular pores (Raoof and

51

Hassanizadeh [14], Joekar-Niasar and Hassanizadeh [4]). The background porosity was assigned a small

52

pore diameter dmin = 0.007 mm. This value is considerably smaller than the pores resolved by CT (0.04

53

mm) and was chosen such that the unsaturated conductivity attains realistic values.

55 56 57

The pore size distribution and the pore topology was adapted to a loess soil which has already been studied by Vogel [15]. We generate a ’structured PNM’ taking account of the measured topology and a ’random PNM’ where the same pores with the same size distribution are connected randomly. The structured

AN US

54

CR IP T

49

and the random grid both have a size of 16 × 16 × 32 nodes. These dimensions are large enough so that the

58

simulations converged to stable results while they are yet small enough to be numerically tractable. The pore

59

size distribution is represented by ten discrete pore size classes. In Table 1 the pore diameters and measured

60

Euler Numbers for the two different PNMs are listed. The definition of the specific Euler number χV referred

61

to as EPC in Table 1 is given by:

62

N +C −H (1) V It is based on the fundamental topological properties of a volume V, which are the number of isolated

63

objects N, the number of redundant connections or loops C often referred to as connectivity or genus and the

64

number of completely enclosed cavities H corresponding to hollow spheres (for more details see Vogel and

65

Roth [12]).

ED

M

χV =

68

(results not shown). This corroborates the well known fact that pore topology is critical for infiltration. We

69

doubled the coordination number from 2.5 (structured PNM) to 5.0 for the random PNM to increase the pore

70

connectivity for the random PNM. As a consequence, the grid constant (i.e. the distance between adjacent

71

nodes) was increased from 0.258 mm (structured PNM) to 0.364 mm and the vertical height of the network

72

from 5.66 mm to 7.98 mm to represent the same porosity for the given number of nodes.

AC

CE

PT

67

When using the same pore size distribution and the same coordination number as for the structured PNM, first simulations indicated that infiltration into the random PNM was blocked due to entrapped gas phase

66

73

74

75

2.2. Calculation of flow To calculate water flow within the pore network as a consequence of some external gradient in water

pressure pw we use Kirchhoffs’ law to guarantee the required mass balance at each node i:

3

ACCEPTED MANUSCRIPT

12 X

qi j = 0

(2)

j=1

77 78

where qi j is the volumetric flow in the pore connecting the nodes i and j, in the following referred to as pore i j. The pressure in the water phase pw is related to the capillary potential pc as follows: pc = pw − pa

79 80

defined as the sum of capillary potential pc and gravitational potential pg :

AN US

The gravitational potential is given by

pg = ρw gh

82 83

(3)

The driving force for vertical infiltration into the porous network is the total hydraulic potential ψw

ψw = pc + pg 81

CR IP T

76

(4)

(5)

with ρw [ML−3 ] the density of water, g [LT −2 ] the gravitational acceleration and h [L] the height above the lower boundary of the network.

84

During transient conditions we need to distinguish three different configurations of phase distribution within an individual pore to calculate qi j : i) completely water saturated, ii) completely gas saturated or iii)

87

partially filled with both phases (i.e. the pore is in the state of infiltration).

88

For water-filled pores (see Fig. 1a)) the flow is obtained by:

ED

85

M

86

qi j = Gi j (ψw,i − ψw, j ),

where the hydraulic conductivity Gi j is given by:

PT

89

(6)

Gi j =

πri4j 8η

(7)

with ri j the radius of the pore i j and η the viscosity. We assume complete wettability of the solid surfaces.

91

During infiltration gas-filled pores are gradually invaded by water. Water can invade a pore with radius ri j from node i (see Fig. 1b)) if and only if all of the following requirements are met:

AC

92

CE

90

93

95

ρw gλi j 1. pc,i + α √ > − 2σ ri j 2 2. There is no water entering from the opposite side of the pore.

96

3. There is a continuous air path from the neighboring node j to the upper or lower boundary of the model.

94

4

ACCEPTED MANUSCRIPT

97

4. No water has entered the pore and ψw,i ≥ ψw, j or the pore is already partly filled.

98

101 102 103

As described in condition 1, the capillary pressure at node i needs to exceed the entry pressure of a given pore described by the Young-Laplace equation. Once the pore is already partially water-filled with a length λi j [L] an additional gravitational component for the water column inside the pore needs to be considered. The weight factor α accounts for the direction of filling, α = {1, 0, −1} for downward, horizontal or upward flow respectively.

CR IP T

99 100

During infiltration of a pore the local flow is described by:

105

! ρw gλi j 2σ qi j = Gmin (ψw,i − ψw, j ) + Gi j pc,i + α √ + (8) ri j 2 The first term represents water flow along the walls of gas-filled pores. The corresponding conductivity

106

Gmin is set to the same value as the background porosity. Such a wall flow is a reasonable assumption

107

and, moreover, it renders the numerical calculation of the pressure field more stable. For the second term

108

describing the pore filling only the capillary potential and not the total hydraulic potential is taken into

109

account because the driving force of the filling process is independent from the gravitational potential at

110

node i. The condition under which a pore is to be filled is the same for all the pores irrespective of their

111

height.

113

For gas-filled and partially filled pores in which bulk flow is not taking place and for ’background porosity’ (see Fig. 1c)) we assume a small water-filled radius rmin so that

M

112

AN US

104

qi j = Gmin (ψw,i − ψw, j )

ED

115

with

Gmin =

4 πrmin 8η

(10)

After infiltration of a given pore, the continuity of the air phase to the lower and upper boundary is

PT

114

(9)

116

updated for each node in the network so that condition 3 can be easily checked in the following time steps.

117

2.3. Initial and boundary conditions Initially the capillary potential is set to a fixed value within the entire network which mimics stationary

119

gravity flow at this potential. The initial water saturation corresponds to the volume of pores that are water

120

filled at this potential. To simulate infiltration we impose a higher but fixed capillary potential at the inlet

121

nodes of the upper boundary. At the lower boundary a gravitational pressure gradient is applied. This is es-

122

tablished by setting the capillary potential at the nodes of the lower boundary equal to the values of the nodes

123

in the second lowest layer. However, this procedure leads to artifacts in case that the hydraulic conductivity

124

of the lower layer (i.e. the bonds between the lowest and the second lowest nodes) is higher compared to

125

the layer above. In this case, the reduced water potential in the second lowest level, which was set equal to

AC

CE

118

5

ACCEPTED MANUSCRIPT

126

the lowest level, induces a feedback so that gravity flow conditions (i.e. a constant matric potential along the

127

vertical axis) are not reached. We solved this problem by changing the pore size distribution of the lower

128

level such that it is smaller than that of the layer above. In this way, the artificial impact of the lower bound-

130

ary on the dynamics within the network is minimized to become negligible. The lateral faces are periodic in the horizontal plane to prevent undesired boundary effects.

a)

b)

CR IP T

129

qi2

qi1

i

i

qij

AN US

qi3

λij

qij

c)

M

i

qij

ED

i

j

CE

PT

qij

AC

Figure 1: Rule based filling mechanisms: a) example for mechanism of water filled pores with bulk flow possible b) example for mechanism of the gradually filling of a gas-filled pore c) examples for mechanism of gas-filled and partly filled pores in which bulk flow is not possible and there exist only flow along the ’background porosity’. Left: pressure level at node i insufficient to fill adjacent bigger pores, right: no bulk flow in pore i j due to entrapped air 131

2.4. Dynamics of infiltration

132

Infiltration was initiated by setting the matric potential at the uppermost nodes, pc,inlet to some value

133

above the current state of the network. In a first step all bonds connected to the upper boundary nodes where

134

ri, j ≥ 2σ/pc,inlet are filled with water. In the following, the dynamics of moving water/air interfaces were 6

ACCEPTED MANUSCRIPT

135

calculated by updating the pressure field and then using the rule based filling conditions after each time-step

136

to calculate fluxes and update water saturation at each pore. The dynamic time-steps are defined by the

137

actual time required to fill the next pore completely. This time-step is calculated with respect to all pores i j

138

containing a moving interface: (

140

For the entire network the flow balance system can be formulated in matrix notation as: A(p)p = b(p)

141

(11)

CR IP T

139

) Va,i j tmin = min qi j Here Va,i j is the air-filled volume of a pore i j and qi j is the flow of the wetting phase into it.

(12)

where A(p) is the n × n matrix containing the conductivities between n individual nodes Gi j and Gmin .

Note that air-filled pores have a low fixed conductivity Gmin . The vector p contains the unknown pressures

143

at the nodes, b is the vector containing the given pressures at the upper and lower boundaries pinlet and poutlet

144

multiplied by the conductances of the bonds connected to these boundaries and the capillary pressures and

145

conductances of pores being filled. A biconjugate gradient method was used to solve the resulting system of

146

linear equations ( Press et al. [16] ). Following the suggestion of Mogensen et al. [10] the calculation of the

147

pressure field was done iteratively by minimizing the function

1 T p Ap − b(p)T p (13) 2 for each time-step. Since f (p) opens upwards and is convex it has an extremum for exactly one vector p

149

M

f (p) =

148

AN US

142

which is the global minimum of f (p). This is obtained based on an iterative approach of the form

150 151

ED

A(pk )pk+1 = b(pk ), k = 0, 1, ...

(14)

When the initial vector p0 is chosen to match the pressure solution of the previous time-step, the number of steps required to minimize f (p) was found to be very small. Solving the linear system of equations a stable and fast solver is necessary. In contrast to the method

153

suggested by Mogensen and Stenby [9] (conjugate gradient method preconditioned by symmetric succes-

154

sive overrelaxation) we used the so called biconjugate gradient method. A mathematical description of the

155

method as well as the computer code of the solver is provided by Press et al. [16]. A set of solvers belonging

156

to the ITPack software package (FORTRAN code - i.e. Jacobi Conjugate Gradient (JCG), Sucessive Over-

157

relaxation (SOR), Symmetric SOR Conjugate Gradient (SSORCG), Reduced System Semi-Iteration (RSSI)

158

solver) was also tested and gave very similar results.

159

3. Continuum scale modeling

AC

CE

PT

152

160 161

To model infiltration at the continuum scale we used the classical Richards equation in its one-dimensional

form 7

ACCEPTED MANUSCRIPT

∂t θ(h) − ∂z K(h)[∂z (h + z)] = 0 162

(15)

where z denotes the vertical coordinate pointing upwards. This equation relies on two hydraulic material

163

properties, the water retention curve θ(h) and the unsaturated hydraulic conductivity function K(h). The

164

volumetric water content is given by θ [-], K [LT−1 ] denotes the hydraulic conductivity and h [L] is the soil

165

matric potential in the form of hydraulic head. The latter is related to capillary pressure by pc = ρw gh. To allow for a direct and consistent comparison between pore scale and the continuum scale simulations

167

we obtained both hydraulic functions from simulations of the PNM. This was done for a stepwise wetting

168

of the PNM to describe the imbibition branch of the hysteretic functions. According to the ten pore size

CR IP T

166

classes of the PNM, this led to 10 discrete points for both functions. For a parametric description of θ(h) we

170

fitted a van Genuchten model (van Genuchten [17]) in its bimodal form according to Durner [18] to increase

171

flexibility. For the hydraulic conductivity function we fitted a spline to the data points. Note that for the

172

scope of this paper it is important to get a good description of the data points simulated by PNM while the

173

number of parameters to do so is not critical.

AN US

169

Finally, based on the parameterized hydraulic functions, Richards equation was solved for the same

175

pressure steps as for the PNM using the software package µPhi (Ippisch et al. [19]). This code uses a cell-

176

centered finite volume scheme for spatial discretization and the solver is based on an implicit Euler scheme

177

in time and a full upwind scheme in space while the time step is adapted automatically.

178

4. Results and discussion

179

4.1. Effective properties

M

174

For both the structured and random PNM we simulated infiltration in 10 pressure steps. After each step

181

we waited for the system to reach hydrodynamic equilibrium (i.e. no changes in local pressures anymore)

182

before initiating the next increase in pressure. We started from a drained model with a pressure of -74 hPa

183

throughout the model slightly lower than the capillary pressure of the smallest measured pore size class.

184

The drainage function which was used to drain the model, was already tested and validated (Vogel and Roth

185

[12]). The first infiltration cycle was simulated by setting the pressures at the top plane of the model (upper

186

boundary) to a pressure of -44 hPa slightly higher than the capillary pressure of the second smallest pore

187

size. During the second step the pressure at the upper boundary was raised to a level of -31 hPa a little bit

188

higher than the capillary pressure of the third smallest pore size and so on.

CE

PT

ED

180

Combined plots for the simulated hydraulic functions of the structured grid in contrast to the random

190

grid are shown in Figs. 2 and 3. The first plot for both models contains the data points of the water content

AC

189

192

over the pressure head obtained by the PNM together with the water retention curve (i.e. capillary pressure saturation relation) as obtained after fitting the van Genuchten parameters of a bimodal model to the output

193

data of the PNM. The fitting function of the bimodal model reads:

191

8

ACCEPTED MANUSCRIPT

Sw =

X

ωi (1 + (αi pc )ni )−(1+1/ni )

(16)

i

Sw = and

195

X

θ − θr θ s − θr

(17)

CR IP T

with the water saturation

194

ωi = 1

(18)

i

with S w : water saturation, θ s : water content at saturation, θr : residual water content, α1 , α2 , n1 ,n2 : van

196 197

Genuchten paramters for submaterial 1 and 2, ω: fraction of the submaterial 1 of the total material.

In the second plot the fitted spline interpolation is shown together with the simulated data point for

199

hydraulic conductivity. We observe significantly higher hydraulic conductivity values for the random PNM. 0.36 0.355 0.35

AN US

198

structured PNM: datapoints fit random PNM: datapoints fit

M

0.34 0.335 0.33

ED

water content θ [-]

0.345

0.325

0.315

-70

CE

0.31 -80

PT

0.32

-60

-50 -40 -30 matric potential pc [hPa]

-20

-10

0

AC

Figure 2: Water retention characteristics simulated by the network Model (symbols) and fitted parameterization (lines) for the structured (blue) and the random PNM (red). The porosity in both networks was fixed to 0.5027

9

ACCEPTED MANUSCRIPT

0

-1

CR IP T

log (permeability K) [cm/h]

-0.5

structured PNM: datapoints fit random PNM: datapoints fit

-1.5

-2.5 0.32

0.325

0.33

AN US

-2

0.335 0.34 0.345 water content θ [-]

0.35

0.355

0.36

Figure 3: Hydraulic conductivities simulated by the network Model (symbols) and spline interpolation (lines) for the structured (blue) and the random PNM (red).

200

There is a considerable difference in hydraulic properties between the random and the structure PNM although the pore size distributions are identical. This is explained by the different topology of the models

202

leading to a different amount of entrapped air within the various pore size classes. This is illustrated in Fig. 4.

203

In the structured grid much more air is entrapped in small pores (water filled at pc < 31 hPa) but less in the

204

larger pores. This is due to the connectivity of the pores as indicated by the Euler numbers (Tab. 1) with

205

lower values for higher connectivity. The connectivity of small pores in the random PNM is increased due

206

to the increased coordination number as explained above while the connectivity of larger pores is slightly

207

increased in the structured PNM reflecting the connectivity of the natural pore space. This explains the

208

different shape of the water retention characteristics.

PT

ED

M

201

The assumption of immiscible fluids leads to a large volume of air entrapments (80-90%) in the larger

210

pores (water filled above pc ≥ 31 hPa). Thus, hydraulic conductivity is mainly governed by the two smaller

CE

209

212

pore size classes which are better connected with an increased water saturation within the random PNM. This is why the random PNM exhibits considerably higher conductivities.

213

We are aware that the enormous volume of air entrapment is not quite realistic for the original soil

214

where pore size and topology was measured. This is rather a consequence of assuming incompressible and

215

immiscible fluids and no additional processes such as dissolution of the gas phase and leaking of bubbles

216

due to local pressure gradients allowing for fluids to escape once they are trapped. However for the scope of

217

this work this is not critical.

AC

211

10

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 4: Entrapped air at different capillary potentials. 218

4.2. Infiltration into the pore network models

To compare the results of pore scale with continuum scale infiltration we simulated infiltration into the

220

two smallest pore size classes of the PNM. Most of the water infiltrates at these pressure steps due to air

221

entrapments in the larger pore size classes. The infiltration process into the entire network from top to the

222

bottom is illustrated in Figs. 5 and 6. Qualitatively this confirms the validity of the rule based pore filling

223

procedure. The infiltration front in the random PNM is remarkably sharper as compared to the structured

224

grid. The reason for this is a more uniform velocity profile due to the random distribution of pores at higher

225

coordination number in contrast to the more clustered topology of the structured PNM.

ED

the behavior of the pore scale properties (i.e. pore size distribution and topology).

AC

CE

227

For a reliable interpretation of the results it is crucial that the PNMs are large enough to actually reflect

PT

226

M

219

11

CR IP T

ACCEPTED MANUSCRIPT

ED

M

AN US

Figure 5: The invasive infiltration process with increasing time (left to right) for the structured grid. First pressure step allowing the two smallest pore size classes to be filled consecutively from top to bottom. Air: red, water: blue.

Figure 6: The same illustration as in previous figure for the random grid.

In Fig. 7 pressure profiles during infiltration are shown. The pressure was calculated as an average over

229

all nodes within a horizontal plane. The curves are smoother as compared to random PNM (shown in Fig. 8)

230

which is due to the sharper infiltration front in the random PNM. Both models converge toward gravity flow

231

and dynamic equilibrium at the end of the infiltration process.

AC

CE

PT

228

12

ACCEPTED MANUSCRIPT

0

-1

CR IP T

depth z in mm

-2

-3

AN US

-4

-5

-6 -75

t=0s t=0.1s t=0.2s t=0.3s t=0.4s t=0.5s t=0.6s t=0.7s t=0.8s t=0.9s t=1.0s t=1.2s t=1.4s t=1.6s t=1.8s t=2.0s t=2.5s t=3.0s t=3.5s t=4.0s t=5.0s t=6.0s t=7.0s t=8.0s t=10s t=20s t=30s t=600s t=4000s

-70

-65

-60 -55 matric potential pc in hPa

-50

-45

-40

AC

CE

PT

ED

M

Figure 7: Pressure profile for the structured grid. The two vertical red lines indicate the capillary pressure of the smallest (d = 0.04mm) and second smallest (d = 0.066mm) pore sizes, at which the pores of these sizes can be filled.

13

ACCEPTED MANUSCRIPT

0

t=0s t=0.1s t=0.2s t=0.3s t=0.4s t=0.5s t=0.6s t=0.7s t=0.8s t=0.9s t=1.0s t=1.2s t=1.4s t=1.6s t=1.8s t=2.0s t=2.5s t=3.0s t=3.5s t=4.0s t=4.5s t=5.0s t=7.0s t=122s

-1

-2

CR IP T

depth z in mm

-3

-4

-5

AN US

-6

-7

-8 -75

-70

-65

-60 -55 matric potential pc in hPa

-50

-45

-40

M

Figure 8: The same plot as in previous figure for the random grid.

For a reliable interpretation of the results it is crucial that the PNMs are large enough so that the observed

233

behavior can actually be related to the pore scale attributes (i.e. pore size distribution and topology). To

234

demonstrate in how far this can be done, we simulated 5 different realizations for both, the structured and

235

the random network. The small standard deviation of the pressure profiles of the different realizations clearly suggests that the volume of the networks is large enough to draw some general conclusions (see Fig. 9).

AC

CE

PT

236

ED

232

14

ACCEPTED MANUSCRIPT

0

t=0s treal. 1=0.1s treal. 2=0.1s treal. 3=0.1s treal. 4=0.1s treal. 5=0.1s treal. 1=1.0s treal. 2=1.0s treal. 3=1.0s treal. 4=1.0s treal. 5=1.0s treal. 1=2.0s treal. 2=2.0s treal. 3=2.0s treal. 4=2.0s treal. 5=2.0s

-1

CR IP T

depth z in mm

-2

-3

-4

-6 -75

-70

-65

AN US

-5

-60 -55 matric potential pc in hPa

0

-2

ED

-4

CE

-6

-40

t=0s treal. 1=0.1s treal. 2=0.1s treal. 3=0.1s treal. 4=0.1s treal. 5=0.1s treal. 1=1.0s treal. 2=1.0s treal. 3=1.0s treal. 4=1.0s treal. 5=1.0s treal. 1=2.0s treal. 2=2.0s treal. 3=2.0s treal. 4=2.0s treal. 5=2.0s

PT

depth z in mm

-3

-5

-45

M

-1

-50

-7

AC

-8 -75

-70

-65

-60 -55 matric potential pc in hPa

-50

-45

-40

Figure 9: Comparisons of the pressure profiles of 5 realizations of the structured (above) and random (below) grid plotted for three different times.

15

ACCEPTED MANUSCRIPT

237

The dynamics of the water content during the same infiltration process is shown in Figs. 10 and 11. It

238

was calculated as an average over the bonds connected to all nodes within a horizontal plane. As for the

239

pressure the water front infiltrates from top toward the bottom. However, the water content profiles show two

241

subsequent waves. First, the smaller pore size class is infiltrated and only when these pores are filled within the entire network the larger pore size class starts to be invaded from the top. This is due to the discrete

240

pore size classes having a fixed pore size each. At the bottom, the water content within the lowest layer is increased because the mean pore size in this layer was reduced to avoid boundary effects as explained in

244

section 2.3..

CR IP T

242 243

245

The infiltration dynamics is similar for the structured and the random PNM while the absolute water con-

246

tent in the random PNM is somewhat higher because of the increased coordination number of this network.

AN US

0

-1

-3

M

depth z in mm

-2

ED

-4

0.325

0.33

0.335 0.34 0.345 water content (volume fraction)

0.35

0.355

Figure 10: Water content profile for the structured grid.

AC

CE

-6 0.32

PT

-5

16

0.36

t=0s t=0.1s t=0.2s t=0.3s t=0.4s t=0.5s t=0.6s t=0.7s t=0.8s t=0.9s t=1.0s t=1.2s t=1.4s t=1.6s t=1.8s t=2.0s t=2.5s t=3.0s t=3.5s t=4.0s t=5.0s t=6.0s t=7.0s t=8.0s t=10s t=20s t=30s t=600s t=4000s

ACCEPTED MANUSCRIPT

0

-1

CR IP T

-2

depth z in mm

-3

-4

-5

AN US

-6

-7

-8 0.32

t=0s t=0.1s t=0.2s t=0.3s t=0.4s t=0.5s t=0.6s t=0.7s t=0.8s t=0.9s t=1.0s t=1.2s t=1.4s t=1.6s t=1.8s t=2.0s t=2.5s t=3.0s t=3.5s t=4.0s t=4.5s t=5.0s t=7.0s t=10.0s t=15.0s t=20.0s t=30.0s t=50.0s t=80.0s t=122s

0.325

0.33

0.335 0.34 0.345 water content (volume fraction)

0.35

0.355

0.36

247

M

Figure 11: Water content profile for the random grid.

4.3. Comparison with the continuum model

Starting from the same initial conditions, the same pressure steps as applied for the PNM simulation

249

were implemented as boundary conditions for 1D Richards equation. We used the water retention curve and

ED

248

the hydraulic conductivity function as shown in Figs. 2 and 3, respectively. To compare the results of the 3D PNM model with those of 1D Richards equation we averaged the water

252

potential and the water content in each horizontal layer of the PNM. The results for the structured grid are

253

shown in Fig. 12. Overall the infiltration profiles are similar, however, there are mainly two remarkable

254

differences. First, the infiltration front within the PNM model is sharper and second, the propagation of the

255

pressure front is initially slower in the PNM and accelerates after the tip of the front reached the bottom.

256

In contrast, the results of Richards equation are more gradual and smoother in space and time. The random

257

PNM (Fig. 13) shows the same behavior but even more pronounced. This is explained by the two discrete

258

pore size classes that are invaded in the PNMs. Each pore size class has a fixed pore diameter so that the

259

critical pressure to infiltrate the larger pores is only reached after the smaller pores were saturated.

AC

CE

PT

251

250

17

ACCEPTED MANUSCRIPT

0

t = 0s t = 1s PNM t = 1s Muphi t = 2s PNM t = 2s Muphi t = 3s PNM t = 3s Muphi t = 4s PNM t = 4s Muphi t = 5s PNM t = 5s Muphi t = 6s PNM t = 6s Muphi t = 7s PNM t = 7s Muphi t = 8s PNM t = 8s Muphi t = 10s PNM t = 10s Muphi t = 20s PNM t = 20s Muphi t = 30s PNM t = 30s Muphi t = 600s PNM t = 600s Muphi t = 4000s PNM t = 4000s Muphi

-1

CR IP T

depth z in mm

-2

-3

AN US

-4

-5

-6 -75

-70

-65 -60 -55 -50 matric potential pc in hPa

-45

-40

AC

CE

PT

ED

M

Figure 12: Pressure profile comparison with Richards equation for the structured grid.

18

ACCEPTED MANUSCRIPT

0

t = 0s t = 1s PNM t = 1s Muphi t = 2s PNM t = 2s Muphi t = 3s PNM t = 3s Muphi t = 4s PNM t = 4s Muphi t = 5s PNM t = 5s Muphi t = 7s PNM t = 7s Muphi t = 122s PNM t = 122s Muphi

-1

-3

CR IP T

depth z in mm

-2

-4 -5

AN US

-6 -7 -8 -75

-70

-65 -60 -55 -50 matric potential pc in hPa

-45

-40

M

Figure 13: Pressure profile comparison with Richards equation for the random grid.

Comparing the water content profiles in Figs. 14 and 15 clearly demonstrates the consecutive infiltration

261

into the two different pore size classes. Individual pores can only be filled in the PNM if the matric potential

262

at the node reaches the capillary pressure of these pores. Obviously this critical capillary pressure is only

263

reached for the larger pore size class when the smaller pores are already infiltrated at the bottom. This

264

leads to the observed stepwise increase in water content and a sharper invasive infiltration front as already

265

discussed for the pressure profiles. Since the water retention characteristic used for the continuum scale

266

model was obtained from PNM-simulations the water content averaged over the entire domain agrees well

267

for the two model concepts. Due to the local heterogeneity within different horizontal planes of the PNM

268

the water content profiles observed in the PNM are less smooth.

AC

CE

PT

ED

260

19

ACCEPTED MANUSCRIPT

0

t = 0s t = 1s PNM t = 1s Muphi t = 2s PNM t = 2s Muphi t = 3s PNM t = 3s Muphi t = 4s PNM t = 4s Muphi t = 5s PNM t = 5s Muphi t = 6s PNM t = 6s Muphi t = 7s PNM t = 7s Muphi t = 8s PNM t = 8s Muphi t = 10s PNM t = 10s Muphi t = 20s PNM t = 20s Muphi t = 30s PNM t = 30s Muphi t = 600s PNM t = 600s Muphi t = 4000s PNM t = 4000s Muphi

-1

CR IP T

depth z in mm

-2

-3

AN US

-4

-5

-6 0.32

0.325

0.33

0.335 0.34 0.345 water content [-]

0.35

0.355

0.36

AC

CE

PT

ED

M

Figure 14: Water content profile comparison with Richards equation for the structured grid.

20

ACCEPTED MANUSCRIPT

0 -1

-3

CR IP T

depth z in mm

-2

t = 0s t = 1s PNM t = 1s Muphi t = 2s PNM t = 2s Muphi t = 3s PNM t = 3s Muphi t = 4s PNM t = 4s Muphi t = 5s PNM t = 5s Muphi t = 7s PNM t = 7s Muphi t = 10s PNM t = 10s Muphi t = 15s PNM t = 15s Muphi t = 20s PNM t = 20s Muphi t = 30s PNM t = 30s Muphi t = 50s PNM t = 50s Muphi t = 80s PNM t = 80s Muphi t = 122s PNM t = 122s Muphi

-4 -5

AN US

-6 -7 -8 0.32

0.325

0.33

0.335 0.34 water content [-]

0.345

0.35

0.355

5. Conclusions

ED

269

M

Figure 15: Water content profile comparison with Richards equation for the random grid.

We simulated water infiltration into a 3-dimensional pore network model representing the statistical

271

properties of the pore size distribution and the pore connectivity of an undisturbed soil. We compared the

272

obtained dynamics of capillary potential and water content to the results of a 1D continuum scale model

273

based on Richards equation. For macroscopically homogeneous media the phenomenology of water dy-

274

namics can be well described based on Richards equation (i.e. Buckingham-Darcy-law). The scope of our

275

pore scale simulation was to explore in how far a discrete pore scale model can reproduce the observed

276

phenomenology at the continuum scale. We came to the following conclusions: • The simulation of infiltration into a 3D pore network by recalculating the entire pressure field after each event of filling single pores leads to infiltration dynamics which are consistent with the classical

AC

278

CE

277

PT

270

279

280

continuum scale approach based on Richards equation.

• The assumption of immiscible and incompressible fluids leads to unrealistic high volumes of en-

281

trapped air during infiltration. In reality there are processes as dissolution of the gas phase and the

282

formation of bubbles to overcome pore blockage by water menisci. It is not obvious how to represent 21

ACCEPTED MANUSCRIPT

283

such process in a physical consistent way within the simplified geometry of a network model. This is

284

expected to be also true for non-cylindrical cross sections of the pores allowing for corner flow.

285

286

6. Acknowledgment We thank DFG for funding within the priority program SPP 1315 and Olaf Ippisch, Kurt Roth for valu-

289

References

291

292 293

294 295

[1] I. Fatt, The network model of porous media. I. Capillary pressure characteristics, Trans. AIME 207 (1956) 144–159.

[2] I. Fatt, The network model of porous media. II. Dynamic properties of a single size tube network, Trans. AIME 207 (1956) 160–163.

AN US

290

CR IP T

288

able support. We are grateful for the assistance of our colleagues Ji-Youn Arns and Samuel Kwame Kumahor.

287

[3] I. Fatt, The network model of porous media. III. Dynamic properties of networks with tube radius distribution, Trans. AIME 207 (1956) 164–181.

[4] V. Joekar-Niasar, S. M. Hassanizadeh, Analysis of fundamentals of two-phase flow in porous media

297

using dynamic pore-network models: A review, Critical Reviews in Environmental Science and Tech-

298

nology 42 (18) (2012) 1895–1976.

300

[5] G. N. Constantinides, A. C. Payatakes, Network Simulation of Steady-State Two-Phase Flow in Consolidated Porous Media, AIChE Journal 42 (1996) 369–382.

ED

299

M

296

[6] B. Ahrenholz, J. T¨olke, P. Lehmann, A. Peters, A. Kaestner, M. Krafczyk, W. Durner, Prediction of

302

capillary hysteresis in a porous material using lattice-Boltzmann methods and comparison to exper-

303

imental data and a morphological pore network model, Advances in Water Resources 31 (9) (2008)

304

1151–1173.

306

307

saturation-interfacial area relationship for porous media, Adv. Water Resour. 32 (2009) 1632–1640.

[8] X. Yang, C. Liu, J. Shang, Y. Fang, V. L. Bailey, A unified multiscale model for pore-scaleflow simulations in soils, Soil Science Society of America Journal 78 (1) (2014) 108–118.

AC

308

[7] M. Porter, M. G. Schaap, D. Wildenschild, Lattice-Boltzmann simulations of the capillary pressure-

CE

305

PT

301

309 310

[9] K. Mogensen, E. H. Stenby, A Dynamic Two-Phase Pore-Scale Model of Imbibition, Transport in Porous Media 32 (1998) 299–237.

311

[10] K. Mogensen, E. Stenby, S. Banerjee, V. A. Barker, Comparison of Iterative Methods for Computing

312

the Pressure Field in a Dynamic Network Model, Transport in Porous Media 37 (1999) 277–301. 22

ACCEPTED MANUSCRIPT

314

315 316

317 318

319 320

321 322

[11] G. Tørå, P.-E. Øren, A. Hansen, A Dynamic Network Model for Two-Phase Flow in Porous Media, Transport in Porous Media 92 (2012) 145–164. [12] H.-J. Vogel, K. Roth, Quantitative morphology and network representation of soil pore structure, Advances in Water Resources 24 (2001) 233–242. [13] A. Raoof, S. M. Hassanizadeh, A New Method for Generating Pore-Network Models of Porous Media, Transport in Porous Media 81 (2010) 391–407, doi:10.1007/s11242-009-9412-3.

CR IP T

313

[14] A. Raoof, S. M. Hassanizadeh, A new formulation for pore-network modeling of two-phase flow, Water Resources Research 48 (2012) 1, doi:10.1029/2010WR010180.

[15] H.-J. Vogel, A numerical experiment on pore size, pore connectivity, water retention, permeability, and solute transport using network models, European Journal of Soil Science 51 (2000) 99–105.

[16] W. H. Press, S. A. Saul A. Teukolsky, W. T. Vetterling, B. P. Flannery, Numerical Recipes in C, The Art

324

of Scientific Computing, chap. Solution of Linear Algebraic Equations, Cambridge University Press,

325

second edn., 83–89, 1992.

327

328 329

[17] M. T. van Genuchten, A closed-form equation for predicting the hydraulic conductivity of unsaturated soils, Soil Sci. Soc. Am. J. 44 (1980) 892–898.

[18] W. Durner, Hydraulic conductivity estimation for soils with heterogeneous pore structure, Water Resour. Res. 30 (1994) 211–223, doi:10.1029/93WR02676.

M

326

AN US

323

[19] O. Ippisch, H.-J. Vogel, P. Bastian, Validity limits for the van Genuchten-Mualem model and impli-

331

cations for parameter estimation and numerical simulation, Adv. Water Resour. 29 (2006) 1780–1789,

332

doi:10.1016/j.advwatres.2005.12.011.

AC

CE

PT

ED

330

23