J_ID: WRCR Customer A_ID: WRCR20978 Cadmus Art: WRCR20978 Ed. Ref. No.: 2013WR014693RR Date: 12-June-14

Stage:

Page: 1

PUBLICATIONS Water Resources Research RESEARCH ARTICLE 10.1002/2013WR014693 Key Points:  Novel regularization of proposed conditions by augmented Lagrangian method  Novel analysis of accuracy and illconditioning of regularizations  Novel evaluation of penalty coefficient as functions of permeability and mesh

Correspondence to: C. Callari, [email protected]

Citation: Abati, A., and C. Callari (2014), Finite element formulation of unilateral boundary conditions for unsaturated flow in porous continua, Water Resour. Res., 50, doi:10.1002/2013WR014693. Received 3 SEP 2013 Accepted 2 JUN 2014 Accepted article online 9 JUN 2014

Finite element formulation of unilateral boundary conditions for unsaturated flow in porous continua

1

A. Abati1 and C. Callari2

2

1

Dipartimento di Ingegneria Civile, University of Rome ‘‘Tor Vergata’’, Rome, Italy, 2DiBT, University of Molise, Termoli, Campobasso, Italy

3

Abstract This paper presents the numerical resolution of unilateral boundary conditions able to effec-

5

tively model several problems of unsaturated flow, as those involving rainfall infiltration and seepage faces. Besides the penalty technique, we also consider the novel regularization of these conditions by means of the more effective augmented Lagrangian method. The performance of the so-obtained finite element method is carefully investigated in terms of accuracy and ill-conditioning effects, including comparisons with analytical solutions and a complete identification of the analogies with the problem of frictionless contact. In this way, we provide a priori estimates of optimal and admissible ranges for the penalty coefficient as functions of permeability and spatial discretization. The proposed method and the estimated coefficient ranges are validated in further numerical examples, involving the propagation of a wetting front due to rainfall and the partial saturation of an aged concrete dam. These applications show that the proposed regularizations do not induce any detrimental effect on solution accuracy and on convergence rate of the employed Newton-Raphson method. Hence, the present approach should be preferred to the commonly used iterative switching between the imposed-flow and the imposed-pressure conditions, which often leads to spurious oscillations and convergence failures.

6

4AQ1

7 8 9 10 11 12 13 14 15 16 17

18 19

ABATI AND CALLARI

1. Introduction

20

In seepage problems, the proper conditions for some boundary portions are often not known a priori. Examples of these problems include rainfall infiltration as well as the presence of a so-called ‘‘seepage face,’’ e.g., in wells, slopes, embankments, and tunnels. Usually, in numerical formulations for partially saturated porous media, these two problems are treated separately, by means of ‘‘atmospheric boundaries’’ [Paniconi and Wood, 1993] and ‘‘variable boundary conditions’’ [Huyakorn et al., 1986], respectively. At the end of each iteration, both these old approaches consider the switching between the imposed-flow and the imposedpressure conditions, depending on problem evolution. It is often observed that the resulting computational methods are unstable, failing to converge, especially if a large number of boundary nodes is simultaneously subjected to the switching of condition [Paniconi and Wood, 1993; Paniconi and Putti, 1994; Geo-Slope, 2012].

21

In spite of this unsatisfying performance, the ‘‘adjustment’’ of boundary conditions is still the most frequent approach in numerical studies [Alonso et al., 2003; Borsi et al., 2006; Bergamaschi et al., 2013] as well as in codes commonly used for problems of subsurface flow [DHI-WASY, 2013] and soil mechanics [Geo-Slope, 2012].

31

As a matter of fact, the two aforementioned boundaries can be effectively treated by a unique general formulation of unilateral conditions on flow, as we proposed with the cooperation of Caro Vargas [2005]. However, in contrast with free-surface approaches to porous flow [Oden and Kikuchi, 1980; Borja and Kishnani, 1991; Zheng et al., 2005], the numerical regularization of unilateral constraints is rarely considered in unsaturated flow formulations. The few available examples are focused only on the penalization of the seepage face condition [Aubry and Ozanam, 1988; Truty and Zimmermann, 2006; Schweizer, 2007; Gerard et al., 2008]. Furthermore, it is still missing a study able to address doubts about penalty sensitivity and ill-conditioning effects introduced by these regularizations. Similar doubts can explain the very limited diffusion of this approach, with respect to the more intuitive switching of boundary conditions.

35

In view of these considerations, herein we present a thorough assessment of the performances exhibited by a numerical formulation of general unilateral boundary conditions on flow. Besides the penalty technique,

44

C 2014. American Geophysical Union. All Rights Reserved. V

ID: padmavathym Time: 17:15 I Path: //xinchnasjn/01journals/Wiley/3B2/WRCR/Vol00000/140266/APPFile/JW-WRCR140266

1

22 23 24 25 26 27 28 29 30

32 33 34

36 37 38 39 40 41 42 43

45

J_ID: WRCR Customer A_ID: WRCR20978 Cadmus Art: WRCR20978 Ed. Ref. No.: 2013WR014693RR Date: 12-June-14

Water Resources Research

Stage:

Page: 2

10.1002/2013WR014693

we also consider the novel regularization of these conditions by means of the more effective augmented Lagrangian method. Both these techniques are implemented in a finite element formulation of the unsaturated flow problem and their efficiency is assessed by means of several numerical examples, including comparisons with analytical solutions and with the formally analogous problem of frictionless contact. In this way, we obtain indications for the proper a-priori setting of the penalty parameters.

46

2. Seepage With Unilateral Constraints

51

For simplicity, we consider the porous space as partially saturated by a single fluid phase and we assume infinitesimal deformations of both the fluid and the solid phase. Furthermore, herein we ignore the coupling between seepage and solid skeleton deformation. However, the straightforward incorporation of the formulation proposed herein in a fully coupled framework [Callari and Abati, 2011] has been positively tested by Callari and Abati [2009]; Callari et al. [2010].

52

2.1. Unsaturated Flow Model The volume of the porous solid is identified by a region X  Rndim for the spatial dimension ndim 51; 2; 3 of the problem. We denote by qw the fluid mass flow and by Mw the fluid mass content, that is, the fluid mass per unit volume of the porous solid. Using a standard notation, the balance of fluid mass can then be written as:

57

_ w 52div qw M

in

X

49 50

54 55 56

58 59 60

(1) 63 62 61

(2)

for the dimensionless coefficient krel of relative permeability and the permeability tensor ksat characterizing the saturated seepage. Consistently with the assumption of infinitesimal deformations, a constant fluid density qw is considered in (2), where we have also denoted by g the gravity acceleration. The relative permeability krel is assumed as dependent on the saturation degree sw, defined as the current fluid volume per unit volume of porous space. The saturation degree is in turn considered as a function of capillary pressure, defined as pc :52pw in single-fluid assumptions.

65 64

The uncoupled formulation considered herein is based on an extension of the fluid content equation by Jacob [1940] to the unsaturated case with compressible solid phase, that is,

71

  2 2 _ w 5qw Cw 1 sw b p_ w M Eoe

Cw 5s2w

b2n nsw0 dsw 1 2n jf dpc js

66 67 68 69 70

72

(3)

for the solid-skeleton oedometric modulus Eoe, the tangent storage modulus

74 73

(4)

and the Biot coefficient b512jsk =js . In (4), we have denoted by sw0 the initial value of the saturation degree and by n the porosity, defined as the void volume for unit volume of the porous solid. We have also used the notation jsk, js, jf for the bulk modulus of solid skeleton, solid phase, and interstitial fluid, respectively.

76 75

2.2. Unilateral Boundary Conditions on Flow Standard conditions on fluid pressure and flow are prescribed at disjoint portions @p X and @q X of the boundary @X of the porous solid, that is,

80

^ w on @p X p w 5p

and

^ wn on @q X qwn 5q

C 2014. American Geophysical Union. All Rights Reserved. V

ID: padmavathym Time: 17:15 I Path: //xinchnasjn/01journals/Wiley/3B2/WRCR/Vol00000/140266/APPFile/JW-WRCR140266

77 78 79

81 82

(5)

for the normal flow qwn :5 qw  n, with n being the outward unit normal to @X. We consider a further boundary region @w X, assuming that @p X; @q X; @w X are mutually disjoint portions whose union is @X

ABATI AND CALLARI

48

53

A generalized Darcy law is used to express the unsaturated flow in terms of the pore pressure field pw: qw 52qw krel ksat ðrpw 2qw gÞ

47

84 83 85

2

J_ID: WRCR Customer A_ID: WRCR20978 Cadmus Art: WRCR20978 Ed. Ref. No.: 2013WR014693RR Date: 12-June-14

Water Resources Research

Stage:

Page: 3

10.1002/2013WR014693

 wn  0 is imposed if the (Figure 1). In @w X, an inflow q pore pressure is less than a prescribed maximum value pmax. Otherwise, an outflow kw occurs as pmax is attained. These requirements are expressed by the following complementarity conditions in terms of the pore pressure ‘‘gap’’ gw: 9 gw :5 qw ðpw 2pmax Þ  0 > > =  wn  0 on @w X kw :5 qwn 2q > > ; kw gw 5 0

87 88 89 90 91

(6)

It can be observed that several problems of practical interest can be solved by using unilateral boundary conditions (6). For example, a potential seepage surFigure 1. Bilateral and unilateral boundary conditions in face can be effectively modeled by setting pmax 5 0 the problem of seepage in a porous solid X, for gw and kw defined by (6)1,2.  wn 5 0, as depicted in Figure 2 for the case of and q a rapid reservoir drawdown. In this problem, zero flow and zero pore pressures should be imposed on the two bank surface portions located above and below the intersection C with the phreatic surface, respectively.

94 93 92 95 96 97 98 99 100 101 102

 wn infiltrates Another example is the rainfall infiltration into an unsaturated slope. The entire given rainfall q the ground if negative water pressures characterize this boundary. However, if ponding is unlikely to occur, the surface pore pressures cannot exceed the atmospheric value and part of the rainfall is shed as runoff. Hence, this situation can be effectively modeled by setting pmax 5 0 in boundary conditions (6).

103

In view of bilateral boundary conditions (5) and of definitions (6)1,2 of gw and kw, the fluid mass balance (1) can be written in weak form as:

107

ð

ð ð _ w dpw dX5 qw  rdpw dX2 M

X

X

 wn q 2 dgw dA2 @w X qw

104 105 106

108

^ wn dpw dA q @q X

ð

ð

kw dgw dA q @w X w

(7)

for all admissible pore-pressure variations dpw 2 V p , with V p :5fdpw : X ! R : dpw 50 on @p Xg and dgw 5qw dpw .

110 109

We remark that the unilateral problem considered herein is formally analogous to the problem of frictionless contact between an elastic body and a rigid obstacle, recalled in Appendix A. This analogy is apparent by comparing the complementarity conditions (6) and (A1) below, as well as the weak forms of the balances of fluid mass (7) and of linear momentum (A2).

112

3. Numerical Regularization

116

In this section, besides the penalty technique, we also consider the novel application of the augmented Lagrangian method for the regularization of unilateral conditions (6). Simple physical interpretations are also provided for both these techniques.

117

In the penalty method, unilateral constraints (6) are fulfilled by imposing

120

kw 5kw hgw i

ABATI AND CALLARI

86

111

113 114 115

118 119

(8)

in the limit kw ! 1 and for the Macaulay operator ‘‘hi’’ returning the positive part of (). In the corresponding numerical formulation, the approximate satisfaction of (6) can thus be obtained by setting a large value of the penalty coefficient kw.

122 121

As indicated by definitions (6)1,2 of kw and gw, the penalty coefficient has the dimensions of the permeability ksat appearing in Darcy law (2) divided by a length. Hence, for a violation hgw i of condition (6)1,

125

C 2014. American Geophysical Union. All Rights Reserved. V

ID: padmavathym Time: 17:15 I Path: //xinchnasjn/01journals/Wiley/3B2/WRCR/Vol00000/140266/APPFile/JW-WRCR140266

123 124

126

3

J_ID: WRCR Customer A_ID: WRCR20978 Cadmus Art: WRCR20978 Ed. Ref. No.: 2013WR014693RR Date: 12-June-14

Water Resources Research

Stage:

Page: 4

10.1002/2013WR014693

Figure 2. Schematic illustration of the application of unilateral boundary conditions (6) to model a potential seepage surface.

multiplier is then locally updated penalizing the calculated gap gwa

coefficient kw can be understood as a surface permeability, localized in @w X and giving the volumetric outflow kw = qw due to a unit value of the porepressure jump hgw i=qw across this boundary.

127

Our numerical formulation of the augmented Lagrangian method is based on the Uzawa iterative algorithm [Laursen, 2002; Wriggers, 2006]. Denoting by a the Uzawa iteration counter, the seepage problem (7) can be solved for a fixed kwa value, and the with a coefficient kw, that is,

133

kwa11 5hkwa 1kw gwa i

128 129 130 131 132

135 136 137 138 139 140

(9)

Complementarity conditions (6) are then fulfilled in both the limits a ! 1 and kw ! 1. In this way, a finite value of kw does not necessarily lead to constraint violation, as it is the case for the penalty method.

143 142 141

The Uzawa update can be understood as the imposition of an outflow at the boundary @w X with surface permeability kw, in order to reduce the violation of the maximum pore pressure constraint.

145

4. Finite Element Formulation

147

We present a finite element method for the resolution of the unsaturated seepage problem described in section 2, focusing the attention on the nodal terms accounting for unilateral constraints (6) by means of penalty and augmented Lagrangian methods.

148

In the generic finite element Xe, the pore pressure pw at point x is approximated by isoparametric interpolation of nodal values pe through a shape function matrix Ne ðxÞ. The evaluation of the pore pressure gradient is based on a standard gradient operator Be : 5 rNe , that is,

151

pwe ðxÞ5Ne ðxÞpe

rpwe 5Be pe

qw 5krel ksat quw with quw :52qw ðBe pe 2qw gÞ

149 150

ð

nel

e51 X e

nel

1A

ð

e51 X e

NTe

152 153

(10) 156 155 154

(11)

Interpolations (10) are used to discretize in space the fluid mass balance (7). The so-obtained ordinary differential equation is written in residual algebraic form and integrated in the generic time step ½tn ; tn11  by a backwardEuler scheme. In this way, denoting by Dt5tn11 2tn the time increment, the residual equation reads: rw 5f ext w 2 A

144

146

The finite-element approximation of the generalized Darcy law (2) then reads:

Mw 2Mw;n dX Dt

159 158 157 160 161

(12)

BTe qw dX50

where symbol A denotes the assembly of contributions of all the nel finite elements and index n 1 1 is omitted to simplify the notation. Equation (12) is given in terms of the following vector of external flows ð nel ^ wn dA fext NTe q w 5sw 2 A e51 @ X q e nel

ð

1ew 2 A

e51 @ X w e

ABATI AND CALLARI

134

163 162 164

(13)  wn NTe q

dA

C 2014. American Geophysical Union. All Rights Reserved. V

ID: padmavathym Time: 17:15 I Path: //xinchnasjn/01journals/Wiley/3B2/WRCR/Vol00000/140266/APPFile/JW-WRCR140266

166 165

4

J_ID: WRCR Customer A_ID: WRCR20978 Cadmus Art: WRCR20978 Ed. Ref. No.: 2013WR014693RR Date: 12-June-14

Water Resources Research

Stage:

Page: 5

10.1002/2013WR014693

where we have used the notations @q Xe :5 @Xe \ @q X and @w Xe :5 @Xe \ @w X. In (13), vector sw includes the given concentrated sources as well as the contributions due to imposed pressures on @p , and nel

ð

ew 52 A

e51 @ X w e

NTe kw dA

167 168 169 170 171 172 173

(14)

is the vector accounting for flow contributions due to unilateral constraints.

175 174

A Newton-Raphson iterative procedure is used to solve algebraic equation (12), whose strong nonlinearity is due to terms accounting for unilateral constraints (14), fluid content equation (3), and generalized Darcy law (11). Thus, at iteration (k 1 1), ðk11Þ ðkÞ the increments Dp5pn11 2pn11 of unknown nodal pressures at tn11 can be calculated by means of the following linearized form of equation (12) obtained using (10 and 11):

178

176 177

179 180 181 182 183 184 185 186 187 188 189

 ðkÞ S rðkÞ 5 Dp (15) 1H2G1P w w Dt

for the standard global matrices of storage S and permeability H, that is, nel

ð

S5 A

e51 X e

NTe

@Mw Ne dX @pw

191 190 192 193

(16) 195 194

nel

ð

H5 A

e51 X e

qw BTe krel ksat Be dX

(17)

respectively, and the operator:

nel

G5 A

ð

e51 X e

BTe

197 196

@krel ksat quw Ne dX (18) @pw

resulting from permeability linearization. To compute the storage Figure 3. Configurations and 1-D finite-element discretizations of the problems: (a) matrix (16), we employ a differeninfiltration in a saturated column; (b) frictionless contact of an elastic solid with a rigid tiation of fluid mass content obstacle. which is consistent with the backward-Euler time integration of (3) used to obtain residual (12) [Lehmann and Ackerer, 1998]. For this time integration, we also consider the mixed form suggested by Celia et al. [1990] to improve mass-balance accuracy. Details of

ABATI AND CALLARI

C 2014. American Geophysical Union. All Rights Reserved. V

ID: padmavathym Time: 17:15 I Path: //xinchnasjn/01journals/Wiley/3B2/WRCR/Vol00000/140266/APPFile/JW-WRCR140266

5

199 198 200 201 202 203 204 205 206

J_ID: WRCR Customer A_ID: WRCR20978 Cadmus Art: WRCR20978 Ed. Ref. No.: 2013WR014693RR Date: 12-June-14

Water Resources Research

Stage:

Page: 6

10.1002/2013WR014693

these developments can be found in the extension to the hydromechanically coupled case presented by Callari and Abati [2009].

207

The further matrix appearing in (15), namely

211

nel

ð

Pw 5 A

e51 @ X w e

NTe

@kw Ne dA @pw

208 209 210

212

(19) 214 213

is due to linearization of the extra term (14) accounting for unilateral constraints. A possible form for ew and Pw can be obtained by means of the penalty method, considering expression (8) for outflow kw: nel

ð

ew 52kw A

e51 @ X w e

215 216 217 218 219 220

NTe hgw i dA (20)

where we have assumed a constant penalty coefficient kw on @w X. So, also in view of the pore-pressure gap definition (6)1, the operator (19) reads: nel

Pw 5kw A

ð

e51 @ X w e

222 221 223 224 225 226

qw Hðgw ÞNTe Ne dA (21)

Figure 4. Convergence profiles of the Uzawa algorithm according to (27) versus the corresponding evolution of constraint violation hgaa i (a: infiltration, b: frictionless contact). Results obtained for different values of the penalty parameter appearing in (9).

for the standard Heaviside function ‘‘H().’’

228 227

As an alternative, the augmented Lagrangian method can be employed to calculate ew and Pw . In this case, the solution in the generic time step requires two nested iterative schemes. In a given Uzawa iteration a, the nonlinear equation (12) is solved by iterative resolution of (15) for a fixed value of outflow kwa . In view of gap definition (6)1, at the generic Newton-Raphson iteration k, the outflow vector (14) and the

230

229

operator (19) then read:

232 233 234 235 236 237 238 239 240 241 242

nel

eðkÞ wa11 52 A

ð

e51 @ X w e

nel

PðkÞ wa11 5kw A

ð

e51 @ X w e

NTe hkwa 1kw gðkÞ wa i dA

(22) 244 243

T qw Hðkwa 1kw gðkÞ wa ÞNe Ne dA

(23)

respectively.

ABATI AND CALLARI

231

246 245

C 2014. American Geophysical Union. All Rights Reserved. V

ID: padmavathym Time: 17:15 I Path: //xinchnasjn/01journals/Wiley/3B2/WRCR/Vol00000/140266/APPFile/JW-WRCR140266

6

J_ID: WRCR Customer A_ID: WRCR20978 Cadmus Art: WRCR20978 Ed. Ref. No.: 2013WR014693RR Date: 12-June-14

Water Resources Research

Page: 7

10.1002/2013WR014693

5. Representative Numerical Simulations

247

For the problem at hand, the performances of both penalty and augmented Lagrangian methods are firstly investigated using a onedimensional element based on a linear interpolation of pore pressures and on a single Gauss quadrature point (section 5.1).

249

Then, we present two-dimensional simulations employing a threenoded triangle with linear interpolation of pore pressures and three Gauss points (sections 5.2 and 5.3). In this case, the penalty regularization of unilateral boundary conditions is implemented in onedimensional linear elements, which are connected to the boundaries @w Xe of the discretized domain. To avoid oscillatory pore-pressure profiles, a two-node Newton-Cotes/ Lobatto quadrature rule is used to calculate (20 and 21) in these onedimensional elements. Both these finite elements have been implemented in code FEAP [Taylor, 2003].

257

For density and unit weight of water, in all the simulations, we assume qw 5103 kg=m3 and cw 510 kN=m3 , respectively. An isotropic permeability is considered, that is ksat 5ksat 1 in (2), for an assumed constant value of saturated permeability ksat. In developments that follow, ksat is given in terms of the corresponding value of hydraulic permeability ksath :5 cw ksat .

275

Furthermore, to avoid oscillatory solutions for very small time increments, a lumped form of storage matrix (16) is used in the numerical simulations presented in this paper.

287

5.1. Infiltration in a Saturated Porous Column In this section, we compare penalty method with augmented Lagrangian procedure in solving a problem of one-dimensional infiltration in a saturated porous continuum. In this comparison, we also consider the formally coincident problem of contact between a linear elastic solid and a rigid obstacle, modeled by the established formulation recalled in Appendix A. The analytical solutions of both these problems are used to estimate the accuracy of flow and stress computed at unilaterally constrained boundaries.

292

Figure 5. Performance of penalty method as a function of penalty parameter (a: infiltration, b: frictionless contact). Runoff/reaction accuracy and constraint violation are assessed by means of error (28) and hga i, respectively. Ill-conditioning is investigated in terms of condition number Nca and minimum tolerance in convergence criterion (31).

ABATI AND CALLARI

Stage:

C 2014. American Geophysical Union. All Rights Reserved. V

ID: padmavathym Time: 17:15 I Path: //xinchnasjn/01journals/Wiley/3B2/WRCR/Vol00000/140266/APPFile/JW-WRCR140266

248

7

250 251 252 253 254 255 256

258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274

276 277 278 279 280 281 282 283 284 285 286

288 289 290 291

293 294 295 296 297

J_ID: WRCR Customer A_ID: WRCR20978 Cadmus Art: WRCR20978 Ed. Ref. No.: 2013WR014693RR Date: 12-June-14

Water Resources Research

Stage:

Page: 8

10.1002/2013WR014693

As illustrated in Figure 3a, a potential  wn is supplied at the water inflow q top of a porous column of height h 5 1.00 m, contained in a box with pervious base, and impervious lateral walls of height h 1 hp 5 1.10 m. As a consequence, the maximum pressure attainable at the column top is pmax 5cw hp , corresponding to the maximum ponding depth hp.

298

In view of the assumed full saturation and of the resulting constant permeability (krel 5 1), the unique nonlinear term in finite-element equation (12) is ew , defined by (14) and due to the unilateral constraint (6) applied to the top of the column  wn and for the aforementioned q pmax. We also assume incompressibility of solid skeleton, solid phase, and water, thus reducing to an elliptic form the problem presented in section 2.1. In this way, besides operator (18), also the storage matrix (16) vanishes in finiteelement equation (15).

308

For the considered permeability value ksath 51026 m=s, the magnitude of the prescribed water inflow  wn 522:1  1023 kg=ðm2 sÞ exceeds q the absolute value of limit inflow qwn 52qw ksath ðh1hp Þ=h 5 21:1 1023 kg=ðm2 sÞ corresponding to the aforementioned maximum pressure pmax. Hence, according to (6)2, the runoff flow

324

kg m2 s (24) is returned to the external environment, thus keeping pw 5pmax at the top of the porous column.

299 300 301 302 303 304 305 306 307

309 310 311 312 313 314 315 316 317 318 319 320 321 322 323

325 326 327 328 329 330 331 332 333

 wn 51:0  1023 kw;ex 5qwn 2q

Figure 6. Comparison between performances of penalty and augmented Lagrangian methods in terms of runoff/reaction errors (28) for different values of penalty parameters (a: infiltration, b: frictionless contact).

It can be observed the coincidence between the particular form of the solving equation (15), with coefficient matrix ðH1Pw Þ, and the equation governing the contact problem depicted in Figure 3b, that is, equation (A6) in Appendix A with coefficient matrix ðK1Pm Þ. In this contact problem, a linear elastic solid of length l 5 2.00 m and Young modulus E 5 200 GPa is constrained at the base. The other end, initially distant gm0 5 21 mm from a rigid obstacle, is loaded by a normal traction t n 5154 MPa, which is greater than the normal stress tn 52Egm0 =l5100 MPa attained in the solid cross sections as a consequence of maximum potential elongation jgm0 j. According to (A1)2, the reaction stress of the rigid obstacle is then: km;ex 5t n 2tn 554 MPa

ABATI AND CALLARI

C 2014. American Geophysical Union. All Rights Reserved. V

ID: padmavathym Time: 17:15 I Path: //xinchnasjn/01journals/Wiley/3B2/WRCR/Vol00000/140266/APPFile/JW-WRCR140266

(25)

8

335 334 336 337 338 339 340 341 342 343 344 345

347 346

J_ID: WRCR Customer A_ID: WRCR20978 Cadmus Art: WRCR20978 Ed. Ref. No.: 2013WR014693RR Date: 12-June-14

Water Resources Research

Stage:

Page: 9

10.1002/2013WR014693

Figure 7. Infiltration in a partially saturated layer: (a) initial pore pressures; (b) boundary conditions; (c) finite element discretization.

For both the infiltration and the contact problems, numerical solutions are obtained through uniform discretizations in 100 one-dimensional finite elements of size dw and dm, respectively (Figure 3). Also for the contact problem, the formulation recalled in Appendix A is implemented in a finite element based on a linear interpolation of displacements and on a single Gauss point.

348

A goal of this work is to estimate the optimal setting for the penalty coefficients of both the augmented Lagrangian and the penalty regularizations. To this purpose, the computations have been repeated for several values of dimensionless ratios  kw =kw and km =km , where the localized values of permeability and stiffness

359

kw :5

ksat dw

 km :5

E dm

350 351 352 353 354 355 356 357 358

360 361 362 363 364 365 366 367 368

(26)

are employed to characterize the considered problem and its spatial discretization. Our choice (26) is suggested by the comparison between the terms of matrices H and Pw used in (15) (or similarly by the comparison between K and Pm ). This setting is also congruent with the physical interpretations of penalty coefficients given in section 3. Furthermore, we note that the use of (26)1 to express the penalty parameter in (8) leads to an outflow expression which is consistent with the regularization of the seepage face condition proposed in the theoretical analysis by Schweizer [2007].

370 369

In the augmented Lagrangian method, the convergence of the Uzawa algorithm at the end of the iteration a is assessed by the criterion:

376

kaa11 2kaa < tolaUZ ka2

for a5w; m

371 372 373 374 375

377

(27)

and for the tolerances tolaUZ . The values of left side of (27) calculated for increasing a are plotted in Figures 4a (a 5 w) and 4b (a 5 m) versus the corresponding constraint violations hgaa i normalized by jga0 j, with gw0 :52qw cw hp and gm0 as defined above. We recall that in a given Uzawa iteration a, the gap gaa is calculated by the Newton-Raphson solution of the global problem (15) for a fixed value of kaa . The updated value kaa11 is then calculated using the Uzawa relation (9). In Figure 4, it can be noted that the correct value of ka

Table 1. Infiltration in a Partially Saturated Layer: Material Parameters

ABATI AND CALLARI

349

Water-saturated permeability

ksath

6  1025

van Genuchten parameter van Genuchten parameter Residual saturation degree Porosity Solid-skeleton volumetric stiffness Solid-skeleton oedometric modulus Solid-phase bulk stiffness Biot coefficient

avg nvg sres w n jsk Eoe js b

0.1 5.6 0.4 0.38 4200 6700 109 1.0

C 2014. American Geophysical Union. All Rights Reserved. V

ID: padmavathym Time: 17:15 I Path: //xinchnasjn/01journals/Wiley/3B2/WRCR/Vol00000/140266/APPFile/JW-WRCR140266

m/s kPa21

kPa kPa kPa

9

379 378 380 381 382 383

J_ID: WRCR Customer A_ID: WRCR20978 Cadmus Art: WRCR20978 Ed. Ref. No.: 2013WR014693RR Date: 12-June-14

Water Resources Research

Stage:

Page: 10

10.1002/2013WR014693

Δ Δ

is approached for decreasing constraint violation hga i, and that the penalty parameter acts as a convergence accelerator for the Uzawa algorithm.

384

In the numerical evaluation of flow and stress at unilaterally constrained boundaries, measures of error can be defined as:

388

Err ðka Þ :5

jka 2ka;ex j ka;ex

385 386 387

389 390 391

for a5w; m (28)

for the exact values kw;ex and km;ex of flow and stress given by the analytical solutions (24) and (25), respectively. In Figures 5a (a 5 w) and 5b (a 5 m), we investigate the performance of the penalty method for a wide range of ratios ka =ka . As expected, the gap violation hga i monotonically decreases for increasing penalty coefficients. The error (28) in runoff/reaction ka decreases until the ratio ka =ka attains optimal values. For larger than optimal penalty coefficients, the observed increase of Err ðka Þ is a consequence of precision loss due to illconditioning effects. This is shown in the same graphs by the increasing condition number, evaluated as Nca : 5jjAa jjp jjA21 a jjp

393 392 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410

for a5w; m (29)

Figure 8. Infiltration in a partially saturated layer: (a) water inflow at top surface for  wn ; (b) pore pressure distributions at different time instants different rainfalls q  wn 520:10 kg=ðm2 sÞ. Results obtained for different temporal and obtained for q spatial discretizations.

in terms of the coefficient matrices Aw 5H1Pw and Am 5K1Pm appearing in equations (15) and (A6), respectively. In (29), ‘‘jj  jjp ’’ denotes the matrix norm corresponding to the p-norm for vectors. In particular, to reduce the computational effort, the p 5 1 norm is considered herein, that is, jjðÞjj1 5 max

1jn

n X

jðÞij j

(30)

i51

Table 2. Partial Saturation of a Concrete Gravity Dam: Material Parameters

ABATI AND CALLARI

Water-saturated permeability

ksath

1028

van Genuchten parameter van Genuchten parameter Residual saturation degree Porosity Solid-skeleton volumetric stiffness Solid-skeleton oedometric modulus Solid-phase bulk stiffness Biot coefficient

avg nvg sres w n jsk Eoe js b

0.05 2.3 0.1 0.12 1:4  107 2:8  107 2:8  107 0.5

C 2014. American Geophysical Union. All Rights Reserved. V

ID: padmavathym Time: 17:15 I Path: //xinchnasjn/01journals/Wiley/3B2/WRCR/Vol00000/140266/APPFile/JW-WRCR140266

m/s kPa21

kPa kPa kPa

10

412 411 413 414 415 416 417 418 419

421 420

J_ID: WRCR Customer A_ID: WRCR20978 Cadmus Art: WRCR20978 Ed. Ref. No.: 2013WR014693RR Date: 12-June-14

Water Resources Research

Stage:

Page: 11

10.1002/2013WR014693

For increasing penalty coefficients, another ill-conditioning effect is illustrated in Figure 5, namely the increasing value of the minimum NR tolerances tolwNR , tolm satisfying the convergence criterion used for the Newton-Raphson procedure, that is,

422 423 424 425 426 427 428

jjraðk11Þ jj < tolaNR jjrð1Þ a jj for a5w; m (31)

As shown in Figure 5, optimal values of penalty coefficients can be chosen in the range 106 4107 ka for both the considered problems. This result is consistent with the recommendation by Taylor [2003] and with the condition obtained by Nour-Omid and Wriggers [1987], taking into account the machine round-off error and the (slight) influence of the problem size. However, we remark that results plotted in Figure 5a suggest that even a wider range of the penalty coefficient can be considered as suitable for flow problems. For example, the range 104 41010 ka would lead to errors in runoff Err ðkw Þ and pressure gap violations hgw i=gw0 less than 1025 .

431 430 429

In Figure 6, we compare the performances of penalty and augmented Lagrangian methods in terms of errors (28) in outflow and stress at unilaterally constrained Figure 9. Partial saturation of a concrete gravity dam: (a) problem geometry and boundary conditions; (b) detail of drainage system with indication of unilateral boundaries. The number of Uzawa constraints. iterations required to attain convergence decreases for increasing values of penalty coefficients. However, for high values of ka , it can be observed a vanishing improvement in runoff/reaction accuracy ensured by the augmented Lagrangian with respect to the penalty method. Hence, for the penalty coefficient used in the augmented Lagrangian method, the range ka 5104103 ka represents a good compromise between runoff/reaction accuracy and convergence rate requirements. It can be noted that this range is consistent with indications reported by Taylor [2003].

450

5.2. Infiltration in a Partially Saturated Layer In this section, we present the application of unilateral boundary conditions (6) to a nonlinear transient problem, namely the rainfall infiltration into an initially nearly dry porous layer. As shown in Figure 7c, this one-dimensional problem is simulated by 236360 triangular finite elements described in section 5. The laws by van Genuchten [1980] are employed to model water retention and relative permeability, that is,

463

8  mvg 1 > > > S ðp Þ5 < w c 11ðavg pc Þnvg h  mvg i2 > > 1=m > : krel ðsw Þ5S1=2 12 12Sw vg w

ABATI AND CALLARI

C 2014. American Geophysical Union. All Rights Reserved. V

ID: padmavathym Time: 17:15 I Path: //xinchnasjn/01journals/Wiley/3B2/WRCR/Vol00000/140266/APPFile/JW-WRCR140266

432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449

451 452 453 454 455 456 457 458 459 460 461 462

464 465 466 467

(32) 469 468

11

J_ID: WRCR Customer A_ID: WRCR20978 Cadmus Art: WRCR20978 Ed. Ref. No.: 2013WR014693RR Date: 12-June-14

Water Resources Research

Stage:

Page: 12

10.1002/2013WR014693

with

470

Sw :5

sw 2sres w and mvg 5121=nvg 12sres w

for pc  0 and the residual saturation degree sres w . The material properties are reported in Table 1, where the values of van Genuchten parameters avg, nvg and sres w are obtained by fitting the experimental data obtained by Topp [1969] for a sandy loam. Simulations presented herein and in the next section employ jf 52  106 kP a for the bulk modulus of water, and an energy criterion to assess convergence of the Newton-Raphson procedure (tolerance value: 10216 ; see [Callari and Abati, 2009] for details).

472 471

The water table is assumed as coincident with the layer base and indifferent to the infiltration process. The initial pore-pressure field is illustrated in Figure 7a. According to (32)1, for the maximum considered value of reference capillary pressure (30 kPa), the residual saturation degree is practically attained. This initial porepressure distribution is almost coincident with that one resulting from a 4 month desaturation period, as shown by simulations of layer desaturation not reported herein for space reasons.

477

The penalty regularization of unilateral constraint (6) is used to simulate the rainfall boundary condition at  wn and no water ponding, that is, pmax 50. the layer top (Figure 7b), for several values of prescribed rainfall q The calculated time evolutions of inflow through the layer top surface are plotted in Figure 8a. The whole  wn infiltrates into the porous solid as far as the pressure at the layer top is below the maximum prerainfall q scribed value (pmax 5 0). As this limit value is attained, the proposed formulation is able to correctly simulate the decrease of the infiltrating flow. As expected, the ponding time required for the attainment of pmax 5 0  wn . decreases for increasing values of rainfall q

482

The propagation of a wetting front in the layer is illustrated in Figure 8b in terms of pore-pressure distributions calculated at different time instants. As time elapses, these profiles get closer to the zero porepressure distribution characterizing the steady state. The infiltration value corresponding to this steady state, that is, jqwn j5qw ksath 50:06 kg=ðm2 sÞ is then the common horizontal asymptote for the curves plotted in Figure 8a.

489

 wn 520:12 kg=ðm2 sÞ, we obtain coincident results for a very fine time As shown in Figure 8a for the rainfall q discretization (10, 620 steps, min 0.05 s, max 4 s, mean 2.6 s) and for a significantly coarser one (326 steps, min 1.6 s, max 128 s, mean 82 s). Also in pore pressure evolutions (not reported herein for space reasons), no apparent loss of accuracy has been observed for the coarser time stepping.

494

As regards the performance of the Newton-Raphson method, the results obtained for the two aforementioned time discretizations show only a slight deterioration of the convergence rate for the larger time steps (from a superquadratic to an asymptotically quadratic rate of convergence). However, for a further coarsening of the time stepping, the Newton-Raphson method exhibits convergence failures. Furthermore, we had to reduce to 1/4 the size of the mentioned coarser time steps to ensure convergence also for the finer mesh (23123120 FEs) obtained by halving the size of the 236360 finite elements. It can be observed that this convergence

498

Figure 10. Partial saturation of a concrete gravity dam: pore pressure distributions at the upstream face of drainage system for different time instants.

ABATI AND CALLARI

(33)

C 2014. American Geophysical Union. All Rights Reserved. V

ID: padmavathym Time: 17:15 I Path: //xinchnasjn/01journals/Wiley/3B2/WRCR/Vol00000/140266/APPFile/JW-WRCR140266

12

473 474 475 476

478 479 480 481

483 484 485 486 487 488

490 491 492 493

495 496 497

499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518

J_ID: WRCR Customer A_ID: WRCR20978 Cadmus Art: WRCR20978 Ed. Ref. No.: 2013WR014693RR Date: 12-June-14

Water Resources Research

Figure 11. Partial saturation of a concrete gravity dam: distributions of pore pressure (a) and saturation degree (b) 103 days after the instantaneous reservoir filling.

ABATI AND CALLARI

Stage:

Page: 13

10.1002/2013WR014693

behavior is consistent with the theoretical results by Radu et al. [2006].

519

A close view of Figure 8a reveals very slight oscillations, which are indifferent to time discretization  wn 52 (see e.g., results for q 0:12 kg=ðm2 sÞ). On the other hand, we have observed that these small inaccuracies completely disappear in solutions obtained with a finer mesh (see  wn 520:10 e.g., the response for q kg=ðm2 sÞ and 23123120 finite elements). However, as illustrated in Figure 8b, the pore pressure profiles obtained with the two considered meshes are practically coincident, thus showing the suitability of the coarser spatial discretization.

522

To assess the performance of the proposed regularization, we have repeated the calculations for a range of the penalty coefficient (kw 5100 4106 m=ðkPasÞ) close to that one evaluated as suitable in section 5.1 (for the problem at hand, it is kw 53:6  1025 m= ðkPasÞ). In the comparison of the corresponding results, no effect of practical relevance has been observed on accuracy and convergence.

540

5.3. Partial Saturation of a Concrete Gravity Dam This section presents the application of the proposed formulation to model the drainage system typically located near the upstream face of concrete gravity dams. As shown in Table 2, we consider the frequent case of an old dam characterized by a relatively high value of the saturated permeability of concrete, due to diffuse damage.

553

Problem geometry and boundary conditions are depicted in Figure 9. To simulate an initially nearly dry dam body, the uniform porepressure distribution pw0 52100

566

C 2014. American Geophysical Union. All Rights Reserved. V

ID: padmavathym Time: 17:15 I Path: //xinchnasjn/01journals/Wiley/3B2/WRCR/Vol00000/140266/APPFile/JW-WRCR140266

13

520 521

523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539

541 542 543 544 545 546 547 548 549 550 551 552

554 555 556 557 558 559 560 561 562 563 564 565

567 568 569 570

J_ID: WRCR Customer A_ID: WRCR20978 Cadmus Art: WRCR20978 Ed. Ref. No.: 2013WR014693RR Date: 12-June-14

Water Resources Research

Page: 14

10.1002/2013WR014693

kPa is imposed for t < 0. According to (32), the corresponding initial value of permeability is krel ksath 56:5  10213 m=s.

571

A rapid reservoir filling is modeled by the instantaneous activation (t 5 0) of a hydrostatic water pressure distribution on the upstream face and of uplift pressures at the dam base. For the uplift distribution, we adopt the typical shape prescribed by dam design codes to account for effects of drains in foundation. In particular, under the drain line, the uplift pressure is set as the 25% of water pressure ph at the dam heel. At  wn 50 and the upstream face of the drainage system, the unilateral constraint (6) is then applied for q

573

( pmax ðzÞ5

ABATI AND CALLARI

Stage:

cw ðhm 2zÞ for z < hm 0

for z  hm

572

574 575 576 577

(34)

where, as illustrated in Figure 9b, we have denoted by z the elevation from the dam base and by hm 50:25 ph =cw the maximum water level in the drains. Consistently with the filling of drains, at the downstream face of the drainage system, a pore pressure distribution pwd(z) coincident with the current upstream profile pwu(z) is imposed below the maximum water level. At higher elevations, the downstream face is set as impervious. As shown in Figure 11, the dam body is discretized by two regular arrangements of 235360 (upstream of drains) and 2310360 (downstream of drains) triangular finite elements. Computations have been successfully performed for a penalty coefficient value (kw 51:0 m=ðkPasÞ) included in the range indicated as admissible in section 5.1 (for the problem at hand, it is kw ’ 1:0  1029 m=ðkPasÞ).

579 578

Pore pressure distributions calculated at upstream face of the drainage line are reported in Figure 10 for different time instants. It is apparent that the maximum pressure profile (34) is progressively attained starting from lower elevations (z < hm) up to the drain portion above the maximum water level hm, with activation of the corresponding outflow (6)2. Fields of pore pressure and saturation degree calculated for t 5 103 days are plotted in Figure 11a and b, respectively. It can be inferred that at this instant, most part of the dam portion located upstream of the drains is almost fully saturated, while the downstream part is practically dry, with the exception of the region close to foundation (Figure 11b). Hence, the proposed finite element treatment of unilateral boundary conditions is able to effectively simulate a behavior feature often observed in aged concrete dams.

587

6. Concluding Remarks

595

The advancements presented in this paper are based on a general formulation of unilateral conditions for flow able to solve problems involving seepage faces and infiltrations. We have developed the novel regularization of these conditions by means of the augmented Lagrangian method. Furthermore, we have assessed the admissibility of the penalty technique for flow problems, providing optimal and admissible ranges for the penalty coefficient as functions of permeability and mesh spacing.

596

The computational method is very effective, as shown by the applications presented herein and in companion papers. The latter publications include 2-D and 3-D simulations of the coupled response of initially saturated banks to reservoir drawdown [Callari and Abati, 2009], the analysis of desaturation induced by dilatancy localization in porous media [Callari et al., 2010] and the 3-D modeling of tunneling [Callari, 2009].

601

These performances fully motivate the selection of properly regularized unilateral constraints as an advantageous alternative to the common switching of boundary conditions, which often induces spurious oscillations and convergence failures.

605

Further developments may include the combination of the proposed formulation with an adaptive time stepping procedure, regulated by a suitable estimation of the truncation error (e.g., see [Kavetski et al., 2001; Clark and Kavetski, 2010] and references therein). This advancement is suggested by the results discussed in section 5.2 and especially by our experience in the large-size 3-D coupled problems aforementioned, where the increase of computational efficiency due to adaptive time stepping would be particularly useful.

608

Appendix A: Signorini Contact Problem

613

We recall a typical finite element formulation of the frictionless contact between a deformable body X and a rigid obstacle. Detailed treatments can be found in [Raous, 1999; Laursen, 2002; Wriggers, 2006]. Standard

614

C 2014. American Geophysical Union. All Rights Reserved. V

ID: padmavathym Time: 17:16 I Path: //xinchnasjn/01journals/Wiley/3B2/WRCR/Vol00000/140266/APPFile/JW-WRCR140266

14

580 581 582 583 584 585 586

588 589 590 591 592 593 594

597 598 599 600

602 603 604

606 607

609 610 611 612

615

J_ID: WRCR Customer A_ID: WRCR20978 Cadmus Art: WRCR20978 Ed. Ref. No.: 2013WR014693RR Date: 12-June-14

Water Resources Research

Stage:

Page: 15

10.1002/2013WR014693

conditions on displacements u and tractions t : 5rn are prescribed at disjoint portions of the boundary @X, that is, u5^ u on @u X and t5^t on @t X, respectively, with n the outward unit normal to @X and r the Cauchy stress tensor, whose normal components are assumed as positive in the tensile case. A further boundary portion @m X can be in frictionless contact with a rigid obstacle, as expressed by the complementarity conditions: 9 gm : 5gm0 1u  n  0 > > = km : 5t n 2tn  0 > > ; km gm 50

@m X

on

s

Denoting by f the body forces and by ‘‘r ðÞ’’ the symmetric part of the gradient tensor, the weak form of linear momentum balance in presence of frictionless contact (A1) reads: ð ð r : rs ðduÞ dX5 f  du dX1 X

ð

X

@m X

620

622 621 623 624 625 626 627

629

(A2)

ð

km dgm dA @m X

ð Xe

BTe r dX50

631 630 632 633 634

(A3)

for the vector of external forces

636 635

nel

fext m 5sm 1 A

e51

ð Xe

ð NTe fdX1

nel

ð

1em 1 A

e51 @ X m e

@t Xe

 NTe ^t dA 1

NTe t n ne

(A4) dA

with @t Xe : 5@Xe \ @t X and ne the outward unit normal to @m Xe : 5@Xe \ @m X. In (A4), vector sm includes the given concentrated forces as well as the contributions due to reactions on @u X, and nel

em 52 A

ð

e51 @ X m e

NTe km ne dA

ðkÞ rðkÞ m 5ðK1Pm Þ Dd

ID: padmavathym Time: 17:16 I Path: //xinchnasjn/01journals/Wiley/3B2/WRCR/Vol00000/140266/APPFile/JW-WRCR140266

639

641 640 642

(A6)

for the nodal displacements d, the standard global matrix of tangent stiffness K, and the operator

C 2014. American Geophysical Union. All Rights Reserved. V

638 637

(A5)

is the vector accounting for contact reactions. A Newton-Raphson iterative procedure may be used to solve equation (A3), whose linearized form reads:

ABATI AND CALLARI

619

628

for all admissible displacement variations du 2 V u , with V u : 5fdu : X ! RndimX : du50 on @u Xg and dgm 5du  n. Using the same notation introduced in section 4, the following standard form of the governing equation, written in residual form, is obtained from substitution of finite element approximations of displacements in balance (A2): nel rm 5fext A m 2 e51

618

^t  du dA @t X

t n dgm dA2

1

617

(A1)

where gm0 and gm are the initial and the current values of the non-positive defined gap between a point on potential contact surface @m X and the rigid body. In (A1)2, we have denoted by tn : 5t  n the normal component of traction vector t on @m X. So, according to (A1)2,3, normal traction tn is equal to the prescribed value t n for a non vanishing gap (gm < 0). If contact is attained (gm 5 0), the corresponding positive-defined reaction stress km is given by (A1)2. As usual, it is assumed that @u X, @t X; @m X are mutually disjoint boundary portions whose union is @X.

ð

616

644 643

15

J_ID: WRCR Customer A_ID: WRCR20978 Cadmus Art: WRCR20978 Ed. Ref. No.: 2013WR014693RR Date: 12-June-14

Water Resources Research

Stage:

Page: 16

10.1002/2013WR014693

nel

ð

Pm 5 A

e51 @ X m e

NTe

  @km ne Ne dA @u

(A7)

due to unilateral term (A5). If the penalization km 5km hgm i is used to regularize (A1), vector (A5) and operator (A7) read: nel

ð

em 52km A

e51 @ X m e

nel

NTe hgm ine dA

Pm 5km A

e51 @ X m e

(A8)

Hðgm ÞNTe ðne

ne ÞNe dA

(A9)

On the other hand, the consideration of the Uzawa update kma11 5hkma 1km gma i leads to the following forms of vector (A5) and operator (A7): nel

ð

nel

ð

e51 @ X m e

652 651 650 653

NTe hkma 1km gðkÞ ma ine dA

(A10)

T Hðkma 1km gðkÞ ma ÞNe ðne ne ÞNe dA

(A11)

e51 @ X m e

PðkÞ ma11 5km A

647

649 648

ð

eðkÞ ma11 52 A

646 645

655 654

658 656 657 Acknowledgments Research supported by MIUR projects on concrete dams (PRIN 2004, code 2004084599-006; PRIN 2007, code 20077ESJAP-002) and on underground storage of CO2 (PRIN 2010-2011, code 2010BFXRHS-004).

ABATI AND CALLARI

References

659

Alonso, E. E., A. Gens, and C. H. Delahaye (2003), Influence of rainfall on the deformation and stability of a slope in overconsolidated clays: A case study, Hydrogeol. J., 11, 174–192. Aubry, D. and O. Ozanam (1988), Free-surface tracking through non-saturated models, in Numerical Methods in Geomechanics, edited by E. Swoboda, pp. 757–763, A. A. Balkema, Innsbruck. Bergamaschi, L., R. Bru, A. Martınez, J. Mas, and M. Putti (2013), Low-rank update of preconditioners for the nonlinear Richards equation, Math. Comput. Model., 57(7–8), 1933–1941. Borsi, I., A. Farina, and M. Primicerio (2006), A rain water infiltration model with unilateral boundary condition: Qualitative analysis and numerical simulations, Math. Method. Appl. Sci., 29, 2047–2077. Borja, R. I., and S. S. Kishnani (1991), On the solution of elliptic free boundary problems via Newton’s method, Comput. Methods Appl. Mech. Eng., 88(2), 341–361. Callari, C. (2009), Three-dimensional effects of ground desaturation due to tunnelling, in Computational Methods in Tunneling—EURO:TUN 2009, edited by G. Meschke et al., pp. 517–524, Aedificatio Publishers. Callari, C., and A. Abati (2009), Finite element methods for unsaturated porous solids and their application to dam engineering problems, Comput. Struct., 87, 485–501. Callari, C., and A. Abati (2011), Hyperelastic multiphase porous media with strain-dependent retention laws, Transp. Porous Media, 86, 155– 176. Callari, C., F. Armero, and A. Abati (2010), Strong discontinuities in partially saturated poroplastic solids, Comput. Methods Appl. Mech. Eng., 199, 1513–1535. Caro Vargas, B. (2005), Formulazione agli elementi finiti di vincoli unilateri nel problema di filtrazione in mezzi parzialmente saturi, MSc thesis, Univ. of Rome ‘‘Tor Vergata’’ and Ecole Nationale des Ponts et Chauss ees, 2005. Celia, M. A., E. T. Bouloutas, and R. L. Zarba (1990), A general mass-conservative numerical solution for the unsaturated flow equation, Water Resour. Res., 26(7), 1483–1496. Clark, M. P., and D. Kavetski (2010), Ancient numerical daemons of conceptual hydrological modeling: 1. Fidelity and efficiency of time stepping schemes, Water Resour. Res., 46, W10510, doi:10.1029/2009WR008894. DHI-WASY Software (2013), FEFLOW 6.2—Finite element Subsurface Flow & Transport Simulation System—User Manual, Berlin. Geo-Slope International Ltd. (2012), Seepage Modeling With SEEP/W, Calgary, Alberta, Canada. Gerard, P., R. Charlier, R. Chambon, and F. Collin (2008), Influence of evaporation and seepage on the convergence of a ventilation cavity, Water Resour. Res., 44, W00C02, doi:10.1029/2007WR006500. Huyakorn, P. S., E. P. Springer, V. Guvanasen, and T. D. Wadsworth (1986), A three-dimensional finite-element model for simulating water flow in variably saturated porous media, Water Resour. Res., 22(13), 1790–1808. Jacob, C. E., (1940), On the flow of water in an elastic artesian aquifer, Trans. AGU, 22, 783–787. Kavetski, D., P. Binning, and S. W. Sloan (2001), Adaptive time stepping and error control in a mass conservative numerical solution of the mixed form of Richards equation, Adv. Water Resour., 24(6), 595–605. Laursen, T. A. (2002), Computational Contact and Impact Mechanics, Springer, Berlin. Lehmann, F., and P. Ackerer (1998), Comparison of iterative methods for improved solutions of the fluid flow equation in partially saturated porous media, Transp. Porous Media, 31, 275–292. Nour-Omid, B., and P. Wriggers (1987), A note on the optimum choice for penalty parameters, Commun. Appl. Numer. Methods, 3, 581–585. Oden, J. T., and N. Kikuchi (1980), Recent advances: Theory of variational inequalities with applications to problems of flow through porous media, Int. J. Eng. Sci., 18, 1173–1284.

660 661 662 663 664 665 666 667 668 669 670 671AQ2 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698

C 2014. American Geophysical Union. All Rights Reserved. V

ID: padmavathym Time: 17:16 I Path: //xinchnasjn/01journals/Wiley/3B2/WRCR/Vol00000/140266/APPFile/JW-WRCR140266

16

J_ID: WRCR Customer A_ID: WRCR20978 Cadmus Art: WRCR20978 Ed. Ref. No.: 2013WR014693RR Date: 12-June-14

Water Resources Research

Stage:

Page: 17

10.1002/2013WR014693

Paniconi, C., and M. Putti (1994), A comparison of Picard and Newton iteration in the numerical solution of multidimensional variably saturated flow problems, Water Resour. Res., 30(12), 3357–3374. Paniconi, C., and E. F. Wood (1993), A detailed model for simulation of catchment scale subsurface hydrologic processes, Water Resour. Res., 29(6), 1601–1620. Radu, F. A., I. S. Pop, and P. Knabner (2006), Newton-type methods for the mixed finite element discretization of some degenerate parabolic equations, in Numerical Mathematics and Advanced Applications, edited by A.B. de Castro et al., pp. 1192–1200, Springer, Berlin. Raous, M. (1999), Quasistatic Signorini problem with Coulomb friction and coupling to adhesion, in CISM Courses and Lectures, No. 384, edited by P. Wriggers and P. Panagiotopoulos, pp. 101–178, Springer, Vienna. Schweizer, B. (2007), Regularization of outflow problems in unsaturated porous media with dry regions, J. Differ. Equations, 237, 278–306. Taylor, R. L. (2003), FEAP—A Finite Element Analysis Program: Version 7.5 Theory Manual, University of California at Berkeley, Berkeley. Topp, G. C. (1969), Soil-water hysteresis measured in a sandy loam and compared with the hysteric domain model, Soil Sci. Soc. Am. Proc., 33(5), 645–651. Truty, A. and T. Zimmermann (2006), Stabilized mixed finite element formulations for materially nonlinear partially saturated two-phase media, Comput. Methods Appl. Mech. Eng., 195, 1517–1546. van Genuchten, M. T. (1980) A closed form equation predicting the hydraulic conductivity of unsaturated soils, Soil Sci. Soc. Am. J., 44, 892– 898. Wriggers, P. (2006), Computational Contact Mechanics, 2nd ed., Springer, Berlin. Zheng, H., D. F. Liu, C. F. Lee, and L. G. Tham (2005), A new formulation of Signorini’s type for seepage problems with free surfaces, Int. J. Numer. Methods Eng., 64(1), 1–16.

699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718

ABATI AND CALLARI

C 2014. American Geophysical Union. All Rights Reserved. V

ID: padmavathym Time: 17:16 I Path: //xinchnasjn/01journals/Wiley/3B2/WRCR/Vol00000/140266/APPFile/JW-WRCR140266

17