A new benchmark semi-analytical solution for density-driven flow in porous media

Accepted Manuscript A new benchmark semi-analytical solution for density-driven flow in porous media Marwan Fahs, Anis Younes, Thierry Alex Mara PII: DOI: Reference:

S0309-1708(14)00080-3 http://dx.doi.org/10.1016/j.advwatres.2014.04.013 ADWR 2194

To appear in:

Advances in Water Resources

Received Date: Revised Date: Accepted Date:

23 January 2014 18 April 2014 20 April 2014

Please cite this article as: Fahs, M., Younes, A., Mara, T.A., A new benchmark semi-analytical solution for densitydriven flow in porous media, Advances in Water Resources (2014), doi: http://dx.doi.org/10.1016/j.advwatres. 2014.04.013

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.

1 2 3 4 5

A new benchmark semi-analytical solution for density-driven flow in

6

porous media

7 8 9 10 11

Marwan Fahs1*, Anis Younes1, Thierry Alex Mara2

12 13 14 15

1

Laboratoire d’Hydrologie de Geochimie de Strasbourg, University of Strasbourg, CNRS, UMR 7517 1 rue Blessig 67000 Strasbourg

16 17

2

Department de Physique, PIMENT, University of La Réunion, Moufia, La Réunion

18 19 20 21 22 23 24 25 26

Submitted to the journal of Advances in Water Resources

27

* Contact person: Marwan Fahs

28

E-mail: [email protected]

29

1

30

Abstract

31

A new benchmark semi-analytical solution is proposed for the verification of density-driven

32

flow codes. The problem deals with a synthetic square porous cavity subject to different salt

33

concentrations at its vertical walls. A steady state semi-analytical solution is investigated

34

using the Fourier-Galerkin method. Contrarily to the standard Henry problem, the cavity

35

benchmark allows high truncation orders in the Fourier series and provides semi-analytical

36

solutions for very small diffusion cases. The problem is also investigated numerically to

37

validate the semi-analytical solution. The obtained results represent a set of new test case high

38

quality data that can be effectively used for benchmarking density-driven flow codes.

39 40 41 42 43 44 45 46 47 48 49 50 51

Key words: density-driven flow, benchmark, semi-analytical solution, Fourier series.

52 53 54 55 56 57 58 59

2

60

1. Introduction

61

Numerical models are considered as irreplaceable tools for the modeling of density-driven

62

flow in porous media [15,22,33,35,39,51,55,57]. They are used for a variety of analyses

63

including seawater intrusion in coastal aquifers, saltwater fingering under sabkha and playa

64

lakes, flow around salt domes as the nuclear waste repositories and saltwater upconing under

65

freshwater lenses. The development and the use of these models require a preliminary

66

validation step to confirm that the nonlinear governing equations are correctly solved. This

67

step is often performed by comparing the results of the numerical models with those of

68

existing benchmark problems. In the literature, a number of benchmarks have been proposed

69

for density-driven flow in porous media [8,18,19,23,24,30,32,34,38,47,52]. The most popular

70

benchmarks [56] are the Henry [23] and the Elder problems [18,19].

71

The Henry problem describes saltwater intrusion into a coastal aquifer. It has been widely

72

used because of the existence of a semi-analytical solution [2,3,4,13,25,36,61]. This problem

73

is considered in [56] as the sole density-driven flow problem with an exact solution. However,

74

the semi-analytical solution of the Henry problem is limited to high values of molecular

75

diffusion coefficient which renders it insensitive to density variations. The worthiness and

76

benefits of this problem for benchmarking density-driven flow codes have been widely

77

studied [16,52,53,56]. All of these studies conclude that this problem is not sufficient to test

78

numerical models.

79

The Elder problem describes the flow induced by a high salinity at the top of a rectangular

80

domain [18,19,27,42,52,53,61] and is considered to be a good benchmark for density-driven

81

codes [56]. It is more sensitive to density variations than the Henry problem. However, due to

82

salt instabilities and fingering phenomena, the Elder problem has no unique solution

83

[50,52,56]. Indeed, contradictory results were reported with either upwelling or downwelling

84

flow in the center of the domain [3,11,21,31,37]. A unique solution to the Elder Benchmark is

85

obtained in [50] using a low Rayleigh number.

86

In this work, we propose a new benchmark semi-analytical solution. The problem deals with a

87

square cavity filled with a saturated porous medium with saline water at one of its vertical

88

walls. It is obtained by recasting the popular thermal porous cavity problem

89

[7,9,10,44,49,54,65] as a variable-density problem where the fluid density is a function of salt

90

concentration. A steady state semi-analytical solution is investigated for both velocity and

91

concentration distributions using the Fourier-Galerkin (FG) method as in the Henry problem

92

[20,23,45,48,66]. To this end, the flow and the transport equations are reformulated in terms 3

93

of stream function and relative concentration. These unknowns are then expanded using

94

appropriate Fourier series truncated at given orders. Finally, applying the Galerkin method

95

with the Fourier terms as trial functions leads to a system of nonlinear algebraic equations

96

having the Fourrier coefficients as unknowns. This system can have a large size for small

97

diffusion coefficient cases that require high truncation orders to achieve stable semi-analytical

98

solutions [45,66]. For the Henry problem, the Fourier expansions induce fairly complex

99

summations that hamper the development of the semi-analytical solution for these cases [14].

100

This difficulty is avoided in the proposed benchmark because all of the terms involving four

101

overlapped summations are reduced to double sums. As a consequence, the semi-analytical

102

solution is investigated for three test cases where the diffusion coefficient is taken to be, ten,

103

one hundred and one thousand times lower than that used in the standard Henry problem.

104

Moreover, an efficient algorithm based on the Powell hybrid method is used to solve the

105

resulting nonlinear system of equations [40,41,43]. The efficiency of this algorithm is

106

improved by an analytical evaluation of the Jacobian matrix.

107

The three test cases are also investigated numerically to assess the worthiness and benefit of

108

the porous cavity problem for benchmarking density-driven flow codes. The numerical

109

simulations are performed using an efficient finite element numerical model developed by

110

Younes and Ackerer [64]. In this model, the flow equation is solved using the Mixed Hybrid

111

Finite Element (MHFE) method [4,58,61,62] and the advection-dispersion transport equation

112

is solved using a combination of the Discontinuous Galerkin (DG) finite element method

113

[46,63] and the Multi-Point Flux Approximation (MPFA) method [1,59-61].

114

The paper is organized as follows: Section 2 is devoted to the description of the benchmark

115

problem and its governing equations. Section 3 is devoted to the development of the semi-

116

analytical solution. In section 4, we show the advantages of the semi-analytical solution for

117

the cavity problem compared to those of the Henry problem. Section 5 briefly describes the

118

numerical solution. In section 6, the semi-analytical and numerical results are presented and

119

discussed in the case of high, small and very small diffusion coefficients. Finally, conclusions

120

are provided in section 7.

121 122

2. Benchmark description and governing equations

123

The benchmark is a synthetic problem inspired from the popular problem of natural

124

convection in porous square cavity [7,9,10,44,49,54,65]. The usual thermal problem is a

125

square box with impermeable boundaries, hot (resp. cold) right (resp. left) vertical wall and

4

126

adiabatic horizontal surfaces. The solute analogous problem, considered in this work, is a

127

square box of size H with no flow boundaries, specified concentration of value one (resp.

128

zero) at the left (resp. right) vertical wall and zero concentration gradient on the horizontal

129

surfaces (Fig. 1). Note that the no flow boundary condition at the left vertical wall is imposed

130

by analogy with the thermal problem. This condition is not physically consistent for high

131

concentrations, since contrarily to the thermal problem, the prescribed boundary concentration

132

constitutes a source of fluid [26].

133

The fluid flow is modeled using the continuity equation and the generalized Darcy’s law.

134

Assuming the Boussinesq approximation to be valid (i.e. the density differences are confined

135

to the buoyancy term), these equations can be written in terms of the equivalent fresh-water

136

head as follows [28]: S

137

q=−

138

∂h + ∇ ⋅q = 0 ∂t

ρ0 g  ρ − ρ0  k  ∇h + ∇z  µ  ρ0 

(1) (2)

139

where S is the mass storage coefficient related to the head changes  L−1  , ρ  ML−3  is the

140

density of the fluid, h is the equivalent freshwater head [ L ] , t is the time [T], q is the

141

Darcy’s velocity  LT −1  , g is the gravity acceleration  LT −2  , µ is the dynamic viscosity of

142

the fluid  ML−1T −1  , k is the permeability  L2  , ρ0  ML−3  is the fresh water density and z

143

is the depth [ L ] .

144

The solute mass transport is governed by the advection dispersion equation:

145

∂C + V ∇C − ∇ ⋅ ( D∇C ) = 0 ∂t

146

where C is the relative solute concentration [ − ] , V = q ε is the fluid velocity, ε is the

147

porosity and D is the dispersion tensor, reduced to molecular diffusion D = Dm I (where I is

148

the identity matrix).

149

The transport equation is coupled to the flow system via the following mixture density

150

equation:

ρ = ρ0 + ( ρ1 − ρ0 ) C

151 152

where ρ1 is the saltwater density.

153

The following initial and boundary conditions are used with the porous cavity problem: 5

(3)

(4)

t = 0 : h = 0, C = 0, 0 ≤ x ≤ H , 0 ≤ z ≤ H t > 0 : qz = 0, ∂C ∂z = 0, z = 0, z = H

154

qx = 0, C = 1, x = 0

(5)

qx = 0, C = 0, x = H 155

where qx and q z are the components of the velocity q in the x and z directions.

156 157 158

3. The semi-analytical solution The steady state continuity equation implies the existence of the stream function defined by:

qx =

159

∂ψ , ∂z

qz = −

∂ψ ∂x

(6)

160

The fresh water head can be eliminated using the curl of Eq. (2). To do so the horizontal and

161

vertical components of this equation are differentiated with respect to z and x respectively,

162

and then subtracted from each other. Using Eq. (6), the resulting equation can be written as

163

follows:  ∂ 2ψ ∂ 2ψ µ +  kg ( ρ1 − ρ 0 )  ∂z 2 ∂x 2

164 165

(7)

Similarly, using Eq. (6), the steady state transport equation simplifies to:  ∂ 2C ∂ 2C  ∂ψ ∂C ∂ψ ∂C Dm ε  2 + 2  = − ∂z  ∂z ∂x ∂x ∂z  ∂x

166 167

 ∂C =  ∂x

(8)

Eqs. (7) and (8) are nondimensionalized using the following variables: Qx =

168

qx H ; Dm

Qz =

ψ qz H x z ; X = ; Z= ; Ψ= Dm H H ε Dm

(9)

169 170

Hence, they can be written as follows:  ∂ 2ψ ∂ 2ψ  ∂C  2 + 2  = Ra ∂Z  ∂x  ∂X

171

 ∂ 2C ∂ 2 C  ∂Ψ ∂C ∂Ψ ∂C −  2 + 2= ∂Z  ∂Z ∂X ∂X ∂Z  ∂X

172

gkH ( ρ1 − ρ0 )

173

where Ra =

174

buoyancy forces and the diffusion effects.

εµ Dm

(10)

(11)

is the Rayleigh number usually defined as the ratio between the

6

175

The resulting system of Eqs. (10) and (11) is solved using the FG method. This method is

176

usually used for accurately solving partial differential equations. The applicability of this

177

method is limited to specific problems with periodic boundary conditions. Indeed, this type of

178

boundary condition is necessary to match the trigonometric function of the Fourier expansion.

179

For the cavity problem, periodic boundary conditions are obtained using the following change

180

of variable:

c = C + X −1

(12)

183

 ∂ 2Ψ ∂ 2Ψ  ∂c − Ra + Ra = 0  2 + 2  ∂X  ∂X  ∂Z

(13)

184

 ∂ 2c ∂ 2c  ∂Ψ ∂c ∂Ψ ∂Ψ ∂c + + =0  2 + 2 − ∂Z  ∂Z ∂X ∂Z ∂X ∂Z  ∂X

(14)

181 182

185

Therefore, Eqs. (7) and (8) become:

This system is subject to the following homogeneous boundary conditions:

186

∂Ψ ∂c = 0, = 0, Z = 0, Z = 1 ∂X ∂Z ∂Ψ = 0, c = 0, X = 0, X = 1 ∂Z

187

These periodic boundary conditions enable the expansion of the stream function Ψ and the

188

concentration c to infinite double Fourier series, truncated at given order as follows:

(15)

Nm Nn

Ψ ( X , Z ) = ∑∑ Am,n sin(mπ Z )sin(nπ X )

189

(16)

m =1 n =1

Nr

Ns

c ( X , Z ) = ∑∑ Br ,s cos( rπ Z )sin( sπ X )

190

(17)

r = 0 s =1

191

where Nm and Nn (resp. Nr and Ns) are the truncation orders for the stream function (resp.

192

concentration) in the x and z directions and Am,n and Br ,s are the Fourier series coefficients

193

for the stream function and the concentration, respectively.

194

Note that the terms sin( mπ Z ) sin( nπ X ) (resp. cos( rπ Z )sin( sπ X ) ) used in the Fourier

195

expansion of the stream function (resp. concentration) identically satisfy the boundary

196

conditions expressed in Eq. (15).

197

Eqs. (16) and (17) are then substituted into Eqs. (13) and (14). The resulting equations are

198

multiplied,

respectively,

by

the

trial

functions

7

ϕ g ,h = 4sin ( gπ Z ) sin ( hπ X )

for

199

( g = 1..Nm, h = 1..Nn )

200

integrated over the square domain. This leads to the following system of nonlinear algebraic

201

equations with Ag , h and Bg , h as unknowns:

202

(

and Θ g ,h = 4 cos ( gπ Z ) sin ( hπ X ) for

)

RFg ,h = π 2 Ag ,h h2 + g 2 −

(

Ra

π

Nr Ns

∑∑ s.B

r ,s

Γ g ,r Γ h,s + Ra

r = 0 s =1

Γ g ,0Γ h,0

π2

( g = 0..Nr , h = 1..Ns ) ,

=0

and

( g = 1..Nm, h = 1..Nn ) (18)

)

RTg ,h = −ε g π 2 h 2 + g 2 Bg ,h + π g Π g ,h 203

−

π2 4

Nm Nn Nr

Ns

∑∑∑∑ Am,n Br ,s ( m.s.λg ,m ,rσ h,n,s + n.r.τ g ,m,rωh,n,s ) = 0 ( g = 0..Nr, h = 1..Ns )

(19)

m =1 n =1 r = 0 s =1

204

where RF and RT are the residuals corresponding to Eqs. (13) and (14) respectively.

205

The coefficients λ , σ , τ , ω and ε g are given by:

206

λg ,m,r = (δ g ,m −r + δ g ,r −m + δ g ,m + r ) ,

σ h ,n ,s = (δ h ,n + s + δ h ,n − s − δ h ,s − n )

(20)

207

τ g ,m,r = (δ g ,r − m + δ g ,m −r − δ g ,m + r ) ,

ω h ,n ,s = ( δ h ,n + s − δ h ,n − s + δ h ,s − n )

(21)

1 → if g ≠ 0 2 → if g = 0

εg = 

208

209

210

211

(22)

The matrices Γ and Π in Eqs. (18) and (19) are: Γi , j =  Ai, j Πi , j =  0

1 − ( −1 )i+ j 1 − ( −1 )i − j + i+ j i− j 1 ≤ i ≤ Nm and 1 ≤ j ≤ Nn

elsewhere

(23)

(24)

212

where δ i , j is the Kronecker delta function.

213

Eqs. (18) and (19) form a system of Nm × Nn + ( Nr + 1) × Ns nonlinear equations. The

214

solution to this system provides the coefficients of the Fourier series and consequently the

215

stream function and the concentration distributions. To solve the nonlinear system of Eqs.

216

(18) and (19), we used an efficient nonlinear solver included in the IMSL library [29]. The

217

solver is based on the modified Powell hybrid algorithm which is an alternative to the

218

Newton’s method [40,41,43]. The algorithm has better rate of convergence than the standard

219

Newton’s method because it estimates the correction using a convex combination of the 8

220

Newton and scaled gradient directions. The algorithm uses the approach of the generalized

221

trust region. Initially, a standard Newton step is performed. If the new solution falls within the

222

trust region, it is used as a trial step in the next iteration. If not, a linear combination of the

223

Newton method and gradient directions is used to predict a new position inside the trust

224

region. The residual is then evaluated at the new position. If the norm of the residual is

225

reduced sufficiently, then the step is accepted and the diameter of the trust region is increased;

226

otherwise, the size of the trust region is decreased and another trial step is computed. The

227

Jacobian matrix is approximated from iteration to iteration using rank one method [12]. With

228

this technique, the entire Jacobian is calculated only at the first iteration or when two

229

successive attempts have failed to reduce the residual. The solver allows for an analytical

230

supply or a numerical (finite difference) calculation of the Jacobian matrix. In this work, the

231

algorithm is supplied with the analytical Jacobian matrix to improve convergence and to

232

reduce computational effort, especially for problems involving a large number of Fourier

233

coefficients. For this purpose, the residuals RF and RT are differentiated with respect to the

234

coefficients Ai , j and Bi , j as follows: ∂RFg ,h

235

∂Ai , j

(

∂RFg ,h

236

∂Bi , j

∂RTg ,h

237

∂Ai, j 238

∂RTg ,h ∂Bi , j

= +π gδ i ,gδ j ,h −

(

)

π2 4

)

= π 2 h2 + g 2 δ g ,iδ h, j

= − Ra

j

π

(25)

Γ g ,i Γh , j

Nr Ns

∑∑ B ( i.s.λ r ,s

g ,i ,r

(26)

σ h , j ,s + j.r.τ g ,i ,rωh, j ,s )

(27)

r = 0 s =1

= −ε gπ 2 h2 + g 2 δ i ,gδ j ,h −

π2 4

Nm Nn

∑∑ A ( m. j.λ m,n

σ h ,n, j + n.i.τ g ,m ,iωh ,n, j )

g ,m,i

(28)

m =1 n =1

239

The size of the resulting nonlinear system (18)-(19) is controlled by the truncation orders of

240

the Fourier series. These orders should be high enough to ensure the stability and accuracy of

241

the solution. In this work, several levels of truncation orders, corresponding to a total number

242

of coefficients ranging from 24 to 90000 coefficients, are investigated for the computation of

243

the solution (see Table 1). The truncation order in the z direction is generally chosen to be

244

greater than in the x direction and the ratio between them is less than 3. This constraint is

245

fixed empirically after several observations to ensure a significant improvement of the semi-

246

analytical solution when going from a level of truncation orders to the next higher level. For

247

all of the investigated truncation levels, the same number of Fourier coefficients is used for

9

248

the concentration and the stream function. This is obtained by imposing Nm=Nr+1 and

249

Nn=Ns.

250 251

4. Advantages of the square porous cavity semi-analytical solution

252

In this section, we showed that contrary to the Henry problem, the semi analytical solution to

253

the cavity benchmark can be applied to cases with small diffusion coefficients. Indeed, the

254

developments shown in the previous section to obtain the semi-analytical solution are similar

255

to those used in the Henry problem [23,45,48,66]. However, the choice of the Fourier

256

expansions terms (Eqs. (16) and (17)) is different because it is dependent on the boundary

257

conditions of the problem. Consequently, the nonlinear system obtained for the Henry

258

problem involves a term with four overlapped summations in the residual of the transport

259

equation (see [45,48,66]). The evaluation of this term requires O ( Nm × Nn × Nr × Ns )

260

operations. This renders the semi-analytical solution impractical for problems requiring high

261

truncation orders as when flow is more affected by buoyancy forces than by diffusion effects.

262

Usually, the ratio between the buoyancy forces and the diffusion effects is measured by the

263

Rayleigh number. The Henry semi-analytical solution has been calculated for Ra ≈ 38 by

264

Henry [23] using 78 Fourier coefficients, by Segol [45] using 138 coefficients and by

265

Simpson and Clement [48] using 203 coefficients. Recently, Zidane et al. [66] succeeded in

266

calculating the semi-analytical solution of the Henry problem for Ra ≈ 175 using 424 Fourier

267

coefficients. The computation of semi-analytical solutions for higher Rayleigh numbers is

268

quickly hampered by the gigantic CPU processing time required for the evaluation of some of

269

the terms in the semi-analytical solution. These difficulties are avoided in the porous cavity

270

benchmark. Indeed, although at first glance similar terms appear in both semi-analytical

271

solutions, a closer look reveals the following simplifications:

272

-

double overlapped summations that require O ( Nr × Ns ) operations (See Appendix A).

273 274

The term involving four overlapped summations can be simplified to four terms of

-

The system of Eqs. (18) and (19) has the particular property that all coefficients Ai , j

275

and Bi , j verifying that i + j is even are equal to zero. This property is demonstrated in

276

Appendix B. This allows for the elimination of half of the unknowns in the final

277

system to be solved.

10

278

-

The analytical Jacobian matrix of the resulting nonlinear system can be calculated

279

without any summation because the sums in Eqs. (27) and (28) reduce to one term

280

when the coefficients λ , σ , τ , ω are substituted into Eqs. (20) and (21).

281 282

This renders the semi-analytical solution for the porous cavity problem more efficient than

283

that for the Henry problem especially for high Rayleigh numbers. Indeed, the previous

284

simplifications make possible the use of a very large number of Fourier coefficients (till

285

90000 ) to calculate high Rayleigh number (till Ra ≈ 38000 ) semi-analytical solutions.

286

5. Numerical solution

287

The semi-analytical solution is validated by comparison with an accurate numerical model

288

developed by Younes and Ackerer [64]. In this model, the coupled flow transport system (1)-

289

(3) is solved on a general triangular mesh using an appropriate spatial discretization for each

290

operator. The flow equation is solved using the MHFE method [4,58,61,62]. This method

291

ensures an accurate and consistent velocity field [17]. It is combined with the mass lumping

292

procedure to avoid unphysical oscillations observed with small time steps [58]. The

293

convection term of the transport equation is solved using the DG method. This method

294

provides accurate solutions with limited numerical diffusion for hyperbolic systems [6,63].

295

The dispersion term of the transport equation is solved using the MPFA method. This method

296

is locally conservative and can treat general irregular grids on anisotropic heterogeneous

297

domains [1,59,60]. In addition, the MPFA and DG methods can be assembled in one system,

298

thus avoiding operator splitting errors [4].

299

The combination of the MFE, DG and MPFA methods was shown to be accurate and robust

300

for modelling density-driven flow problems [4]. In this model the fluid flow and transport

301

equations are solved sequentially using the non-iterative time stepping scheme based on local

302

truncation error control [64].

303

The porous cavity benchmark is simulated using a uniform triangular mesh obtained by

304

subdividing square elements into four equal triangles (by connecting the center of the square

305

to its four nodes). Several levels of grid refinement are used to study the sensitivity of the

306

numerical solution to the grid size. For a given grid level n, the number of squares is

307

( 25 × n )

308

corresponds to a total number of triangles ranging from 2500 (level 1) to 90000 (level 6).

309

Transient simulations are performed until a long non-dimensional time to ensure steady

310

solutions.

2

. Six levels of grid refinement are used with n ranging from 1 to 6, which

11

311

6. Results and discussion

312

In this section, the semi-analytical and numerical solutions of the porous cavity problem are

313

investigated for three test cases with different diffusion coefficients. In these test cases, the

314

diffusion coefficient is taken to be: ten, one hundred and one thousand times lower than that

315

used in the Henry problem. Aside from the diffusion coefficient, all the other parameters are

316

similar to those used in the Henry problem (Table 2).

317

For each test case, the semi-analytical and numerical solutions are calculated using different

318

levels of truncation orders and grid refinement. The accuracy of the semi-analytical and the

319

numerical solutions are assessed regarding several parameters related to the flow and the mass

320

transfer processes. As is customary for cavity problems, the streamlines, the horizontal (resp.

321

vertical) velocity profile at the mid-section X = 0.5 (resp. Z = 0.5 ), the maximum values of

322

the dimensionless horizontal velocity Qxmax at the mid-section X = 0.5 and that of the

323

dimensionless vertical velocity Qzmax at the mid-section Z = 0.5 are used for the assessment of

324

the flow process. For the mass transfer process, the accuracy of the solution is assessed using

325

the principal (0.25, 0.5 and 0.75) concentration contours and the average Sherwood number

326

(analogous of the Nusselt number in heat transfer). This number represents the rate of mass

327

transfert across the enclosure. The local Sherwood number is given by the following

328

relationship [49]: Sh =

329 330

∂C ∂X

(29) X =0

The average Sherwood number is defined as follow: 1

Sh = ∫ Sh.dZ

331

(30)

0

332

The first simulation results show that the average Sherwood number is much more sensitive

333

than concentration contours when varying the level of truncation orders in the semi-analytical

334

solution or when varying the level of grid refinement in the numerical solution. Therefore, the

335

value of the average Sherwood number as well as Qxmax and Qzmax are used to check the

336

convergence of both the semi-analytical and numerical solutions.

337

The relative sensitivity of these parameters to the truncation orders at level (k) is calculated as

338

follows:

339

SenYTr ,k =

1 ∆Y 1 Y k − Y k −1 = Y k ∆Nc Y k Nc k − Nc k −1

12

(31)

340

where Y refers to the parameter of interest ( Sh , Qxmax and Qzmax ), ∆Y and ∆Nc represent

341

the variations of Y and Nc between the truncation levels k and k-1 and Nc is the total

342

number of Fourier coefficients given by Nc = Nm × Nn + Nr × Ns .

343

Similarly, the sensitivity of these parameters to the grid refinement at a certain level (n) is

344

defined as follows: SenYGd , n =

345

1 ∆Y 1 Y n − Y n −1 = Y n ∆Ne Y n Ne n − Nen −1

(32)

346

where ∆Y and ∆Ne represent the variations of Y and Ne between the grid levels n and n-1

347

and Ne is the total number of elements given by Ne = 4 ( 25 × n ) for the grid level n .

348

In all the following test cases, we considered the convergence of the semi-analytical (resp.

349

numerical) solution to be achieved when the relative sensitivity of the parameters Sh , Qxmax

350

and Qzmax to truncation orders (resp. grid refinement) for becomes less than 10 −6 .

351

It should be noted that the accuracy of the proposed semi-analytical solution is assessed not

352

only with global quantities (such as the concentration contours) but also with local parameters

353

(such as the Sherwood number and the maximum velocities). This makes the semi-analytical

354

solution suitable for studying the accuracy of new numerical methods and schemes developed

355

in the context of density-driven flow [3,4,5,13,37,61,64].

356

2

-

Test case 1: high diffusion coefficient ( Ra = 380 )

357

In this test case, the diffusion coefficient is taken to be ten times lower than the one used in

358

the standard Henry problem. The corresponding Rayleigh number is approximately 380. The

359

semi-analytical solution is calculated for 7 levels of truncation orders corresponding to a total

360

number of Fourier coefficients ranging from 24 to 6000 (as in Table 1). The convergence of

361

the nonlinear solver is reached for all levels of truncation orders starting with zero as the

362

initial values for the Fourier coefficients. An analytical evaluation of the Jacobian matrix is

363

used because its finite difference approximation is impractical for a large number of Fourier

364

coefficients. As an example, for level 6 with 3000 coefficients, the numerical evaluation of

365

the Jacobian requires more than 1 day of CPU time, whereas, the same matrix can be

366

calculated analytically in 15 seconds.

367

The semi-analytical solutions calculated for the first three levels show strong unphysical

368

oscillations. These oscillations vanish when using higher truncation orders. The convergence

369

of the semi-analytical solution is achieved at the level 7 involving 6000 Fourier coefficients

370

because all sensitivities are less then 10−6 . The principal concentration contours (0.25, 0.5, 13

371

0.75), the streamlines and the velocity profiles obtained with the truncation levels 4 through 7

372

are coincident. The converged average Sherwood number, maximum velocities ( Qxmax and

373

Qzmax ) and their sensitivity to truncation orders are provided in Table 3.

374

The test case is then simulated using the numerical model for three levels of grid refinement.

375

The results show that the average Sherwood number and the maximum velocities are weakly

376

sensitive to grid size. At grid level 3, all sensitivities are less than 10−6 . Therefore, the

377

convergence of the numerical solution is considered to be achieved at this level. Note that

378

similar concentration contours, streamlines and velocity profiles are obtained with the three

379

grid levels. The converged average Sherwood number, maximum velocities ( Qxmax and Qzmax )

380

and their sensitivity to grid size are depicted in Table 4.

381

The converged numerical and semi-analytical solutions are then compared regarding their

382

maximum velocities, average Sherwood number, concentration contours, streamlines and

383

velocity profiles. Very closely converged maximum velocities and average Sherwood

384

numbers are observed in Tables 3 and 4. Moreover, the results of Figs. 2a, 2b and 3a show a

385

strong agreement between the converged numerical (level 3) and semi-analytical (level 7)

386

concentration contours, streamlines and velocity profiles, respectively. The semi-analytical

387

velocity field is illustrated in Fig. 2a. This figure shows that saltwater invades the domain

388

from the left side to the right side. The density effects generate counter-clockwise circulating

389

flow with relatively slow motion in the central region (Fig. 2b). The flow promotes saltwater

390

intrusion near the bottom of the cavity. Fig. 2a shows that the transition from saltwater to

391

fresh water regions in the cavity is relatively smooth. The thickness of the transition zone,

392

defined by the distance between the 0.1 and the 0.9 concentration contours, increases from

393

0.59 near the bottom (at Z= 0.05) to 0.9 at the center of the domain (at Z = 0.5).

394

To distinguish buoyancy effects from diffusion effects, the simulation of a tracer (without

395

density contrast) is performed using grid level 3. The concentration contours obtained are

396

vertical and the resulting Sherwood number is 1.08 (instead of 3.77 for the coupled problem).

397

This suggests that; in this test case, approximately 71% of the mass transfer rate in the domain

398

is due to density effects.

399

-

Test case 2 small diffusion coefficient ( Ra = 3800 )

400

In this test case, the diffusion coefficient is taken to be one hundred times lower than the one

401

used in the Henry problem. The corresponding Rayleigh number is approximately 3800. In

402

this case, high truncation orders are required to guarantee the stability of the semi-analytical 14

403

solution. This solution is calculated for levels 5 through 11 of the truncation orders, which

404

correspond to a total number of Fourier coefficients ranging from 1000 to 36000 (as in Table

405

1). The resulting systems of equations are strongly nonlinear and special attention should be

406

given to the choice of the initial values to ensure the convergence of the nonlinear solver. To

407

this aim, we proceed in the following manner. For each truncation level, the semi-analytical

408

solution is first evaluated for Dm = 18.86 ×10−6 m2 S −1 . Then, the resulting Fourier coefficients

409

are used as initial guess to find the solution for a lower diffusion coefficient. This procedure is

410

repeated successively by decreasing the diffusion coefficient until reaching the desired value.

411

The semi-analytical results show that the average Sherwood number and the maximum

412

velocities are more sensitive to truncation level than in the first test case. The convergence of

413

the semi-analytical solution is reached at level 11 involving 36000 Fourier coefficients.

414

Similar concentration contours, streamlines and velocity profiles are obtained with truncation

415

levels 8 through 11.The converged average Sherwood number, maximum velocities and their

416

sensitivity to the truncation orders are given in Table 3.

417

It is relevant to recall here that, contrary to the Henry problem, the use of a large number of

418

Fourier coefficients was possible because of the simplification of the term involving four

419

overlapped summations. Indeed, the evaluation of the residual with four overlapped

420

summations in the case of 36000 Fourier coefficients requires 324 × 106 operations. This takes

421

more than 3 days of CPU time. The evaluation of the same term in the simplified form

422

(involving double overlapped summations) requires 18000 operations and takes only 10 sec

423

of CPU time.

424

The test case is then simulated numerically using four levels of grid refinement. The results

425

show that the average Sherwood number and the maximum velocities are more sensitive to

426

grid size than in the previous test case. They also show that, contrary to the previous test case,

427

the concentration contours, streamlines and velocity profiles are slightly sensitive to grid size.

428

The converged numerical solution is obtained at the grid level 4 because all sensitivities

429

become less than 10−6 . This convergence is confirmed using the concentration contours,

430

streamlines and velocity profiles because grid levels 3 and 4 lead to the same results. The

431

converged average Sherwood number, maximum velocities and their sensitivity to grid size

432

are depicted in Table 4.

433

The converged numerical (grid level 4) and semi-analytical (truncation level 11) solutions are

434

then compared regarding their average Sherwood number, maximum velocities, concentration

435

contours, streamlines and velocity profiles. A strong agreement is displayed between the semi 15

436

analytical and numerical results in Tables 3 and 4 and in Figs. 3b, 4a and 4b. The velocity

437

field is illustrated in Fig. 4a. This figure shows that the transition zone from saltwater to

438

freshwater regions becomes narrower than in the previous test case especially on the vertical

439

walls near the top left and bottom right corners. The thickness of the transition zone, increases

440

from 0.2 near the bottom (at Z = 0.05) to 0.97 at the center of the domain (at Z = 0.5).

441

The tracer simulation using grid level 4 results in an average Sherwood number of 1.47

442

(instead of 15.89 for the coupled problem). This means that approximately 91% of the final

443

mass transfer rate in the domain is due to density effects. -

444

Test case 3: very small diffusion coefficient ( Ra = 38000 )

445

In this test case the diffusion coefficient is taken to be one thousand times lower than the one

446

used in the Henry problem. The corresponding Rayleigh number is approximately 38000. The

447

semi-analytical solution is calculated for levels 8 through 14 of the truncation orders, which

448

correspond to a total number of Fourier coefficients ranging from 9000 to 90000 (as in Table

449

1). As with the other test cases, the convergence of the semi analytical solution for each

450

truncation level is reached by starting with the Fourier coefficients obtained for

451

Dm = 18.86 ×10−6 m 2 s −1 and then decreasing successively until reaching the desired value

452

(D

453

The results show that the average Sherwood number and the maximum velocities are very

454

sensitive to truncation order. The convergence of the semi-analytical solution is considered to

455

be achieved at the truncation level 14 involving 90000 Fourier coefficients. The concentration

456

contours, streamlines and velocity profiles for levels 12, 13 and 14 show very strong

457

agreement. Note that the computation of the semi-analytical solution at this level is made

458

possible due to the simplification detailed in Appendix B which allows for the elimination of

459

half of the Fourier coefficients from the final system to be solved.

460

The converged average Sherwood number, the maximum velocities and their sensitivity to the

461

truncation level are provided in Table 3.

462

Numerical simulations are performed using 6 levels of grid refinement. The results show that

463

the Sherwood number, maximum velocities, concentration contours, streamlines and velocity

464

profiles are very sensitive to grid size. The numerical solution at grid level 6 involving 90000

465

triangles is considered to be the converged solution. The grid levels 5 and 6 provide similar

466

concentration contours, streamlines and velocity profiles.

467

The converged numerical (grid level 6) and semi-analytical (truncation level 14) Sherwood

468

numbers and maximum velocities show strong agreement because their difference is less 3%

m

)

= 18.86 ×10−9 m2 s −1 .

16

469

for Sh , 2% for Qxmax and 6% for Qzmax . The converged semi-analytical and numerical

470

concentration contours, streamlines and velocity profile are given in the Figs. 5a, 5b and 3c,

471

respectively. These Figures depict excellent matching between the numerical and the semi-

472

analytical solutions. Fig. 5a illustrates the velocity field. It shows that, in this case, the

473

concentration distribution becomes very steep and the transition zone is very narrow

474

especially when close to the vertical walls near the top left and bottom right corners. In these

475

regions, the thickness of the transition zone is only 0.05 (at Z = 0.05) due to sharp

476

concentration profile. The thickness of the transition zone reaches 0.995 at the center of the

477

domain (at Z = 0.5).

478

The average Sherwood number produced by the tracer simulation with at grid level 6 is 2.47

479

(as opposed to 52.12 for the coupled problem). This means that; in this last case,

480

approximately 95% of the final mass transfer rate across the domain is due to density effects.

481 482

In the Fig. 6a we studied the variation of the truncation error of the FG method as a function

483

of the number of Fourier coefficients for the three studied test cases. This error is defined to

484

be the mean of the relative truncation errors on Sh , Qxmax and Qzmax . At the truncation level (k)

485

it is given by: Er k =

486

1 Y k − Y cv Y cv Y = Sh ,Qxmax ,Qzmax 3

∑

(33)

487

where Y cv represents the converged semi-analytical value of the parameter of interest Y .

488

Similarly, the Fig. 6b shows the variation of the mean relative truncation error of the FE

489

solution as a function of the number of unknowns for the three studied test cases. The Fig. 6a

490

(respectively 6b) allows for estimating the required number of coefficients (respectively

491

unknowns) to obtain a desired accuracy with the FG (respectively FE) method. These figures

492

show that the number of required coefficients (or unknowns) increases when the Rayleigh

493

number increases. The slope of the three curves in the Fig. 6a is more important than that in

494

the Fig. 6b suggesting that the convergence of the FG method is much faster than that of the

495

FE method.

496

Finally, for future inter-codes comparison and/or code validations, we provided the local

497

Sherwood number at X=0 for the three test cases. These profiles are given in Fig. 7.

498 499 500 17

501

7. Summary and conclusion

502

In this work, a new benchmark semi-analytical solution has been developed for the

503

verification of density-driven flow numerical codes. The proposed benchmark addresses a

504

square cavity filled with a porous medium that is subject to different concentrations at its

505

vertical walls.

506

The semi-analytical solution has been developed using the Fourier-Galerkin (FG) method as

507

with the standard Henry problem [23]. To this aim, the governing equations were formulated

508

in terms of relative concentration and stream function which are expanded to a double set of

509

Fourier series truncated at a given order. The resulting highly nonlinear system of algebraic

510

equations was solved using the modified Powell’s algorithm with an analytical evaluation of

511

the Jacobian.

512

Contrary to the Henry problem, the semi analytical solution of the cavity benchmark supports

513

cases with very small diffusion coefficients, which are known to require a large number of

514

Fourier coefficients. Several levels of truncation orders have been investigated for the semi-

515

analytical solution of three test cases where the diffusion coefficient is taken to be, ten, one

516

hundred and one thousand times lower than that used in the standard Henry problem.

517

The three test cases have also been investigated numerically using an efficient numerical

518

model where the flow equation is solved using the MHFE method and the convection-

519

diffusion transport equation is solved using a combination of the DG and MPFA methods.

520

Several levels of grid refinement have been used to check the convergence of the numerical

521

solution.

522

In the first test case, corresponding to a Rayleigh number of 380, the convergence of the semi-

523

analytical solution was reached using 6000 Fourier coefficients. The average Sherwood

524

number, maximum velocities, concentration contours, streamlines and velocity profiles are

525

weakly sensitive to grid size. The convergence of the numerical solution is achieved with

526

22500 elements. The converged numerical and semi-analytical solutions show excellent

527

agreement regarding the average Sherwood number, maximum velocities, concentration

528

contours, streamlines and velocity profiles. For this Rayleigh number, the transition from

529

saltwater to freshwater regions is relatively smooth and the transition zone is relatively wide.

530

In the second test case, corresponding to a Rayleigh number of 3800, high truncation orders

531

are required to guarantee the stability of the semi-analytical solution. In this case, special

532

attention should be given to the initial guess to ensure the convergence of the nonlinear

533

solver. The converged semi-analytical solution is obtained with 36000 Fourier coefficients,

18

534

whereas, the converged numerical solution is achieved with 40000 elements. The average

535

Sherwood number and the maximum velocities become sensitive to grid size contrary to the

536

concentration contours, streamlines and velocity profiles which remain only slightly sensitive.

537

A strong agreement is observed between semi-analytical and numerical solutions. At this level

538

of Rayleigh number, the transition zone becomes narrower than in the first case especially

539

near the vertical walls.

540

In the last test case, corresponding to a Rayleigh number of 38000, the buoyancy effects

541

become dominat. Therefore, the concentration distribution becomes very steep and the

542

transition zone becomes very narrow especially close to the vertical walls near the top left and

543

the bottom right corners. The converged semi-analytical and numerical solutions are achieved

544

with 90000 Fourier coefficients and 90000 triangles, respectively. The difference between the

545

converged numerical and semi-analytical is around 3%, 2% and 6% for the average Sherwood

546

number and the maximum horizontal and vertical velocities, respectively. Concentration

547

contours, streamlines, velocity profiles and average Sherwood number are very sensitive to

548

grid size. Therefore, accurate simulation of this test case may be considered as a real

549

challenge for computational models.

550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 19

568

Appendix A: Development of the nonlinear term

569

The term involving four overlapped summations can be written as the summation of the two

570

following sub-terms Nm Nn Nr Ns

Term1g ,h = ∑∑∑∑ Am,n Br ,s ( m.s.λg ,m,rσ h,n ,s )

571

(A.1)

m =1 n =1 r = 0 s =1 Nm Nn Nr

Ns

Term2 g ,h = ∑∑ ∑∑ Am ,n Br ,s ( n.r.τ g ,m ,r ωh,n ,s )

572

(A.2)

m =1 n =1 r = 0 s =1

573

Using the expressions of the coefficients λ and σ as given in Eq. (20), Eq (A.1) can be

574

written as follows: Nm Nn Nr

Ns

Nm Nn Nr

Ns

Term1g ,h = ∑∑∑∑ m.s.Am,n Br ,sδ g ,m− rδ h,n+ s + ∑∑∑∑ m.s.Am,n Br ,sδ g ,m− rδ h,n− s m =1 n =1 r = 0 s =1 Nm Nn Nr

m =1 n =1 r = 0 s =1

Ns

−∑∑∑∑ m.s.Am,n Br ,sδ g ,m− rδ h ,s −n m =1 n =1 r = 0 s =1 Nm Nn Nr Ns

Nm Nn Nr

Ns

m =1 n =1 r = 0 s =1

m =1 n =1 r = 0 s =1

+∑ ∑∑∑ m.s.Am,n Br ,sδ g ,r − mδ h ,n + s + ∑∑∑∑ m.s.Am,n Br ,sδ g ,r −mδ h,n− s 575

Nm Nn Nr

Ns

(A.3)

−∑∑∑∑ m.s.Am,n Br ,sδ g ,r − mδ h ,s −n m =1 n =1 r = 0 s =1 Nm Nn Nr Ns

Nm Nn Nr Ns

m =1 n =1 r = 0 s =1

m =1 n =1 r = 0 s =1

+∑∑∑∑ m.s.Am,n Br ,sδ g ,m+ rδ h,n+ s + ∑∑∑∑ m.s.Am,n Br ,sδ g ,m+ rδ h,n− s Nm Nn Nr

Ns

−∑∑∑∑ m.s.Am,n Br ,sδ g ,m+ rδ h,s− n m =1 n =1 r = 0 s =1

576

Let us consider the first term in Eq. (A.3). Using the relations δ g ,m −r = δ r ,m − g and δ h ,n + s = δ h −n ,s

577

and the property

∑∑∑∑ A

i,j

i

j

k

Bk ,lδ i ,k δ j ,l = ∑∑ Ai , j Bi , j , this term can be simplified to be:

l

i

Nm Nn Nr Ns

∑∑∑∑ m.s.A

578

m,n

j

Nm Nn

Br ,sδ g ,m−rδ h,n + s = ∑∑ m.( h − n ) .Am,n Bm− g ,h− n

m =1 n =1 r = 0 s =1

(A.4)

m =1 n =1

579

When the same procedure is applied to all terms in Eq (A.3), Term1g ,h can be written as

580

follows: Nm Nn

581

I II III  Term1g ,h = ∑∑ mAm,n ( h − n ) Βm,n,g ,h + ( n − h ) Βm,n ,g ,h − Βm,n ,g ,h ( h + n ) 

(A.5)

m =1 n =1

582

583

I II III where Βm,n,g ,h , Β m,n,g ,h and Β m,n,g ,h are given by:

Β Im,n,g ,h = ( Bm− g ,h− n + Bg + m,h− n + Bg −m,h− n ) 20

(A.6)

584

Β IIm,n,g ,h = ( Bm− g ,n− h + Bg + m,n− h + Bg −m,n− h )

(A.7)

585

Β mIII,n,g ,h = ( Bm − g ,h+ n + Bg + m ,h+ n + Bg − m ,h+ n )

(A.8)

586

The term Term2 g ,h is simplified in the same manner. Hence, it can be written as follow: Nm Nn

Term2 g ,h = ∑∑ nAm,n  ΒmIV,n,g ,h ( g + m ) + ΒVm,n ,g ,h ( m − g ) − ΒVIm,n ,g ,h ( g − m ) 

587

(A.9)

m =1 n =1

588

where ΒmIV,n ,g ,h , ΒVm ,n ,g ,h and ΒVIm,n ,g ,h are given by:

589

Β IV m,n,g ,h = ( Bg + m,h − n − Bg + m ,n − h + Bg + m ,h + n )

(A.6)

590

ΒVm,n,g ,h = ( Bm− g ,h− n − Bm− g ,n−h + Bm− g ,h+ n )

(A.7)

591

ΒVIm,n,g ,h = ( Bg −m ,h− n − Bg −m ,n− h + Bg − m,h + n )

(A.8)

592

Appendix B: Reduction of the nonlinear system dimension

593

The residual of the flow equation ( RFg ,h ) is given by Eq. (18). In this equation the term

594

Γ g ,r Γ h ,s (in the double overlapped summation) can be written as follows:

Γ g ,r Γ h,s = 595

1 − ( −1 )g + r 1 − ( −1 )h + s 1 − ( −1 )g + r 1 − ( −1 )h− s 1 − ( −1 )g −r 1 − ( −1 )h+ s + + g+r h+s g +r h−s g −r h+s

1 − ( −1 )g − r 1 − ( −1 )h − s + g−r h−s

(B.1)

596

By a simple development of the Eq. (B.1), we can show that the term Γ g ,r Γ h ,s vanishes in the

597

case where g + h is even and r + s is odd. In the same way we can show that the last term

598

(Γ

599

( RF ) corresponding to

600

the Fourier coefficient Ag ,h and Bg ,h that verify that g + h is even. The same treatment of the

601

residual of the transport equation

602

coefficients Ag ,h and Bg ,h leading to g + h being even can be eliminated for the final system

603

and solved separately. All these coefficients are equal to zero because the independent

604

corresponding system is homogeneous.

g ,0

Γ h,0 ) in Eq. (18) drops out if g + h is even. This means that the residual equations

g ,h

g + h being even are all homogenous. These equations involve only

( RT ) g ,h

leads to a final conclusion that the Fourier

605 606 21

607

References

608

[1] Aavatsmark I. An introduction to multipoint flux approximations for quadrilateral grids.

609

Computat. Geoscien. 2002;6:404-432. http://dx.doi.org/10.1023/A:1021291114475.

610

[2] Abarca E, Carrera J, Sanchez-Vila X, Dentz M. Anisotropic dispersive Henry problem,

611

Adv. Water Resour. 2007; 30(4):913–926. http://dx.doi.org/10.1016/j.advwatres.2006.08.005.

612

[3] Ackerer P, Younes A, Mosé R. Modeling variable density flow and solute transport in

613

porous medium: 1. Numerical model and verification. Transp. in Porous Media 1999;

614

35(3):345-373. http://dx.doi.org/10.1023/A:1006564309167

615

[4] Ackerer P, Younes A. Efficient approximations for the simulation of density driven flow

616

in porous media. Adv Water Res 2008; 31:15-27.

617

[5] Albets-Chico X, Kassinos S. A consistent velocity approximation for variable-density

618

flow and transport in porous media. Journal of Hydrology 2013; 507:33–51.

619

[6] Arnold DN, Brezzi F, Cockburn B, Marini LD. Unified analysis of discontinuous Galerkin

620

methods for elliptic problems. SIAM J. Numer. Anal. 2002;5:1749–79.

621

[7] Baez E, Nicolas A. 2D natural convection flows in tilted cavities: Porous media and

622

homogeneous fluids Int. J. Heat and Mass Trans. 2006;49:4773–4785.

623

[8] Bakker M, Oude Essink GHP, Langevin CD. The rotating movement of three immiscible

624

fluids – a benchmark problem. J Hydrol 2004;287:270–8.

625

[9] Baytas AC. Entropy generation for natural convection in an inclined porous cavity. Int. J.

626

Heat Mass Transfer 2000;43:2089–2099.

627

[10] Bejan A. On the boundary layer regime in a vertical enclosure filled with a porous

628

medium. Lett. Heat Mass Transfer 1979;6:93–102.

629

[11] Boufadel MC, Suidan MT, Venosa AD. Numerical modelling of water flow below dry

630

salt lakes: effect of capillarity and viscosity. J Hydrol 1999:55–74.

631

[12] Broyden CG. A Class of Methods for Solving Nonlinear Simultaneous Equations.

632

Mathematics of Computation (American Mathematical Society) 1965;19(92):577–593.

633

doi:10.2307/2003941.

634

[13] Buès M, Oltéan C. Numerical simulations for saltwater intrusion by the mixed hybrid

635

finite element method and discontinuous finite element method. Transport Porous Med

636

2000;40(2):171–200.

637

[14] Costa B. Spectral Methods for Partial Differential Equations. Cubo Mathematical Journal

638

2004;6:1-32. 22

639

[15] Diersch H-JG. FEFLOW reference numerical, finite element subsurface flow and

640

transport simulation system. Report in WASY institute for water resource planning and

641

systems research, Berlin, 2002.

642

[16] Diersch H-JG, Kolditz O. Variable-density flow and transport in porous media:

643

approaches and challenges. Adv Water Resour 2002;25:899–944.

644

[17] Durlofsky L. Accuracy of mixed and control volume finite element approximations to

645

Darcy velocity and related quantities. Water Resour Res 1994;30:965–73.

646

[18] Elder JW. Steady free convection in a porous medium heated from below. J Fluid Mech

647

1967;27:29–48.

648

[19] Elder JW. Transient convection in a porous medium. J Fluid Mech 1967;27:609–23.

649

[20] Forbes LK. Surface waves of large amplitude beneath an elastic sheet. Part 2. Galerkin

650

solution, J. Fluid Mechanics 1988;188:491-508.

651

[21] Frolkovic P, De Schepper H. Numerical modelling of convection dominated transport

652

coupled with density driven flow in porous media. Adv Water Res2001;24:63–72.

653

[22] Gossel W, Sefelnasr A, Wycisk P. Modelling of paleo-saltwater intrusion in the northern

654

part of the Nubian Aquifer System, Northeast Africa. Hydrogeol J 2010;18:1447–63.

655

[23] Henry HR. Effects of dispersion on salt encroachment in coastal aquifers. Water-Supply

656

Paper 1613-C. US Geological Survey; 1964.

657

[24] Henry JW, Hilleke JB. Exploration of multiphase fluid flow in a saline aquifer system

658

affected by geothermal heating. Washington DC: Bureau of Engineering Research, Report

659

150–118, US Geological Survey Contract 14–08-0001-12681, NTIS, Publication PB234233;

660

1972.

661

[25] Herbert AW, Jackson CP, Lever DA., Coupled groundwater flow and solute transport

662

with fluid density strongly dependent upon concentration. Water Resour Res 1988; 24:1781–

663

95.

664

[26] Hidalgo J., Carrera J. Effect of dispersion on the onset of convection during CO2

665

sequestration, J. Fluid Mech. 2009;640:441–452

666

[27] Holzbecher EO. Modeling density-driven flow in porous media: principles, numeric,

667

software. Springer, Heidelberg, Germany.

668

[28] Huyakorn PS, Andersen PF, Mercer JW, White HO. Saltwater intrusion in aquifers:

669

Development and testing of a three-dimensional finite element model, Water Resour. Res.

670

1987;23:293-312.

671

[29] IMSL. Math/Library. Houston, TX: Visual Numerics; 1997.

23

672

[30] Johannsen K, Kinzelbach W, Oswald S, Wittum G. The salt pool benchmark problem –

673

numerical simulation of saltwater upcoming in a porous medium. Adv Water Resour

674

2002;25:335–48.

675

[31] Kolditz O, Ratke R, Diersch HJ, Zielke W. Coupled groundwater flow and transport: 1.

676

Verification of variable-density flow and transport models. Adv Water Res 1998;21:27–46.

677

[32] Konikow LF, Sanford WE, Campbell PJ. Constant concentration boundary conditions:

678

Lessons from the HYDROCOIN variable-density groundwater benchmark problem. Water

679

Resour Res 1997;33:2253–61.

680

[33] Langevin CD, Thorne DT, Dausman AM, Sukop MC, Guo W. SEAWAT version 4: a

681

computer program for simulation of multispecies solute and heat transport. US Geo. Sur.

682

techniques and methods, book 6, chap. A22, US Geo Sur., Reston, 2007.

683

[34] Langevin CD, Dausman AM, Sukop MC. Solute and heat transport model of the Henry

684

and Hilleke laboratory experiment. Ground Water 2010;48:757–70.

685

[35] Nishikawa T, Siade AJ, Reichard EG, Ponti DJ, Canales AG, Johnson TA. Stratigraphic

686

controls on seawater intrusion and implications for groundwater management, Dominguez

687

Gap area of Los Angeles, California,USA. Hydrogeol J 2009;17:1699–725.

688

[36] Oldenburg CM, Pruess K. Dispersive transport dynamics in a strongly coupled

689

groundwater-brine flow system. Water Resour Res 1995;31:289–302.

690

[37] Oltean C, Buès MA. Coupled groundwater flow and transport in porous media. A

691

conservative or non-conservative form? Transport Porous Med 2001;44(2):219–46.

692

[38] Oswald SE, Kinzelbach W. Three-dimensional physical benchmark experiments to test

693

variable-density flow models. J Hydrol 2004;290:22–42.

694

[39] Pool M, Carrera J. Dynamics of negative hydraulic barriers to prevent seawater intrusion.

695

Hydrogeol J 2010;18:95–105.

696

[40] Powell MJD. An efficient method for finding the minimum of a function of several

697

variables without calculating derivatives. Comp. J. 1964;7:155–162.

698

[41] Powell MJD. On the calculations of orthogonal vectors, Comp. J. 1968;1:302–304.

699

[42] Prasad A, Simmons CT. Using quantitative indicators to evaluate results from variable-

700

density groundwater flow models. Hydrogeol J 3006;13(5-6):905-914.

701

[43] Press WH, Teukolsky SA, Vetterling WT, Flannery BP. Numerical Recipes: The Art of

702

Scientific Computing, 3rd ed., chap. 9-10, Cambridge University Press., New York, 2007.

703

[44] Saeid NH, Mohamad AA. Periodic free convection from a vertical plate in a saturated

704

porous medium, non-equilibrium model. Int. J. Heat Mass Transf. 2005;48:3855–3863.

24

705

[45] Segol G. Classic Groundwater Simulations Proving and Improving Numerical Models,

706

Prentice-Hall, Old Tappan, N. J, 1994.

707

[46] Siegel P, Mose R, Ackerer P, Jaffré J. Solution of the advection diffusion equation using

708

a combination of discontinuous and mixed finite elements. Int. J. Numer. Methods In Fluids

709

1997;24;595–613.

710

[47] Simmons CT, Narayan KA, Wooding RA. On a test case for density-dependent

711

groundwater flow and solute transport models: The salt lake problem. Water Resour Res

712

1999;35:3607–20.

713

[48] Simpson MJ, Clement TP. Theoretical analysis of the worthiness of the Henry and Elder

714

problems as benchmarks of density-dependent groundwater flow models. Adv. Water Resour.

715

2003;26:17-31.

716

[49] Sivasankaran S, Kandaswamy P,

717

density fluids in a porous cavity. Transp Porous Med 2008;71:133–145.

718

[50] van Reeuwijk M, Mathias SA, Simmons CT, Ward JD. Insights from a pseudo-spectral

719

approach

720

http://dx.doi.org/10.1029/2008WR00742.

721

[51] Voss CI, Provost AM. SUTRA, a model for saturated-unsaturated variable density

722

ground water flow with energy or solute transport. US Geol. Surv. Open-File Rep 02-4231,

723

2002, 250.

724

[52] Voss CI, Simmons CT, Robinson NI. Three-dimensional benchmark for variable-density

725

flow and transport simulation: matching semi-analytic stability modes for steady unstable

726

convection in an inclined porous box. Hydrogeol J 2010;18:5–23.

727

[53] Voss CI, Souza WR. Variable-density flow and solute transport simulation of regional

728

aquifers containing a narrow freshwater–saltwater transition zone. Water Resour Res.

729

1987;23:1851–66.

730

[54] Walker KL, Homsy G.M. Convection in a porous cavity. J. Fluid Mech. 1978;87:449–

731

474.

732

[55] Watson TA, Werner AD, Simmons CT. Transience of seawater intrusion in response to

733

sea

734

http://dx.doi.org/10.1029/2010WR00956.

735

[56] Werner A, Bakker M, Post V, Vandenbohede A, Lu C, Ashtiani B, Simmons C, Barry

736

DA. Seawater intrusion processes, investigation and management: Recent advances and future

737

challenges. Adv. Wat. Resour. 2013; 51:3–26.

to

level

the

rise.

Elder

Ng CO. Double diffusive convection of anomalous

problem.

Water

Water

Resour.

25

Resour

Res.

Res

2010;

2009;45:W04416.

46:

W12533.

738

[57] Yechieli Y, Shalev E, Wollman S, Kiro Y, Kafri U. Response of the Mediterranean and

739

Dead Sea coastal aquifers to sea level variations. Water Resour Res 2010;46:W12550.

740

http://dx.doi.org/10.1029/2009WR00870.

741

[58] Younes A, Ackerer P, Lehmann F. A new mass lumping scheme for the mixed hybrid

742

finite element method. Int. J. Num. Meth. Eng. 2006;67:89-107.

743

[59] Younes A, Fontaine V. Hybrid and Multi Point Formulations of the Lowest Order Mixed

744

Methods for Darcy's Flow on Triangles. Int J Numer Meth in Fluids 2008;58:1041-1062.

745

[60] Younes A, Fontaine V. Efficiency of Mixed Hybrid Finite Element and Multi Point Flux

746

Approximation methods on quadrangular grids and highly anisotropic media. Int. J. Numer.

747

Meth. in Eng. 2008;76(3):314-336.

748

[61] Younes, A, Fahs M, Ahmed S. Solving density flow problems with efficient spatial

749

discretizations and higher-order time integration methods, Advances in Water Resources

750

2009; 32:340-352.

751

[62] Younes A, Ackerer P, Delay F. Mixed finite element for solving diffusion-type

752

equations. Rev. Geophys. 2010;48:doi:10.1029/2008RG000277.

753

[63] Younes A, Fahs M, Ackerer P. An efficient geometric approach to solve the slope

754

limiting problem with the Discontinuous Galerkin method on unstructured triangles, Int. J.

755

Numer. Meth. In Bio. Eng. 2010;12:1824-1835.

756

[64] Younes, A, Ackerer P. Empirical versus time stepping with embedded error control for

757

density-driven

758

doi:10.1029/2009WR008229.

759

[65] Zheng W, Robillard L, Vasseur P. Convection in a square cavity filled with anisotropic

760

porous media saturated with water near 4°C. Int, J. Heat and Mass Trans. 2001;44:3463-3470.

761

[66] Zidane A, Younes A, Huggenberger P, Zechner E. The Henry semi-analytical solution

762

for

763

2012;doi:10.1029/2011WR011157.

saltwater

flow

in

porous

intrusion

with

media,

Water

reduced

764

26

Resour.

Res.

dispersion. Water

2010;46:W08523,

Resour.

Res.

765

Figure captions

766

Fig. 1. Domain for the porous cavity benchmark.

767 768

Fig. 2. Test case 1: (a) Concentration contours and velocity field and (b) streamlines.

769 770

Fig. 3 . Dimensionless velocity profiles for the three Test cases: horizontal velocity (Qx) at

771

mid-section X=0.5 and vertical velocity (Qz) at mid-section Z=0.5.

772

Fig. 4. Test case 2: (a) Concentration contours and velocity field and (b) streamlines.

773 774

Fig. 5. Test case 3: (a) Concentration contours and velocity field and (b) streamlines.

775 776

Fig. 6. Truncation error variation: (a) for the FG method and (b) for the numerical model (Nb

777

is the number of unknowns in the numerical model).

778 779

Fig. 7. Local Sherwood number at X=0

780 781 782 783 784

27

Figure 1

Fig. 1. Domain for the porous cavity benchmark.

1

Figure 2

(a)

(b)

Fig. 2. Test case 1: (a) Concentration contours and velocity field and (b) streamlines.

Figure 3

(a)

20

40

15

Test Case 1 S-Anal (level-Tr7) Num (level-Gd 3)

10

Test Case 1 S-Anal (level-Tr 7) Num (level-Gd 3)

20

5

Qx

Qz

0

0

-5 -20

-10 -15

-40

-20 0,0

0,2

0,4

0,6

0,8

0,0

1,0

0,2

0,4

(b)

0,6

0,8

1,0

0,8

1,0

0,8

1,0

X

Z 600

100 80

40

Test Case 2 S-Anal (level-Tr 11) Num (level-Gd 4)

400

Test case 2 S-Anal (level-Tr 11) Num (level-Gd 4)

60

200

20

Qx

0

Qz

0

-20 -40

-200

-60 -400

-80 -100 0,0

0,2

0,4

0,6

0,8

-600

1,0

0,0

0,2

0,4

Z

0,6

X

400

(c)

6000

Test Case 3 S-Anal (level-Tr 14) Num (level-Gd 6)

300

Test case 3 S-Anal (level-TR 14) Num (level-Gd 6)

200

4000 2000

100

Qx

Qz

0

0

-100

-2000

-200

-4000

-300

-6000

-400 0,0

0,2

0,4

0,6

Z

0,8

1,0

0,0

0,2

0,4

0,6

X

Fig. 3. Dimensionless velocity profiles for the three Test cases: horizontal velocity (Qx) at mid-section X=0.5 and vertical velocity (Qz) at mid-section Z=0.5

Figure 4

(a)

(b)

Fig. 4. Test case 2: (a) Concentration contours and velocity field and (b) streamlines.

1

Figure 5

(a)

(b)

Fig. 5. Test case 3: (a) Concentration contours and velocity field and (b) streamlines.

1

Figure 6

10-1

Er

Er

10-1

10-2

10-2

Test Case 1 Test Case 2 Test Case 3

Test Case 1 Test Case 2 Test Case 3

10-3 103

104

105

10-3 104

105

Nc

(a)

106

Nb

(b)

Fig. 6. Truncation error variation: (a) for the FG method and (b) for the numerical model (Nb is the number of unknowns in the numerical model)

Figure 7

Test Case 1 Test Case 2 Test Case 3

100

Sh

10

1

0,1 0,0

0,2

0,4

0,6

0,8

Z

Fig. 7. Local Sherwood number at X=0

1,0

785 786

Table 1: Truncation levels used for the computation of the semi-analytical solution (Nc is the total number of Fourier coefficients).

Tr-level

Nm=Nr+1

1

24

3

Nn=N s 4

2

60

5

6

3

300

10

15

4

600

15

20

5

1000

20

25

6

3000

25

60

7

6000

40

75

8

9000

50

90

9

16000

80

100

10

24000

100

120

11

36000

120

150

12

60000

150

200

13

74800

170

220

14

90000

180

250

Nc

787 788

28

Table 2: Parameters for the porous cavity problem.

789 790

Square size

H=1m

Permeability

kx = k y = k = 1.0204 ×10−9 m2

Porosity Molecular Diffusion (case 1)

ε = 0.35 Dm = 18.86 ×10−7 m 2 s −1

Molecular Diffusion (case 2)

Dm = 18.86 ×10−8 m 2 s −1

Molecular Diffusion (case 3)

Dm = 18.86 ×10−9 m 2 s −1

Fresh water density

ρ0 = 1000 kg/m3

Salt water density

ρ1 = 1025 kg/m3

Gravity

g = 9.81 m/s2

Viscosity

µ = 10−3 kg.m-1 .s−1

791

29

Table 3: Semi-analytical average Sherwood number, dimensionless maximum velocities and

their sensitivities to the truncation orders for the three test cases.

Test case 1

Level (k) 7

Nc 6000

2

11

3

14

Sh

,k SenTr Sh

Qxmax

,k SenQTrmax x

,k SenQTrmax

48.53

6.1×10-7

z

-3.77

1.0×10

-7

20.90

9.3×10

36000

-15.89

6.5×10-7

89.07

1.9×10-7

567.33

4.3×10-7

90000

-52.12

3.0×10-7

391.30

9.9×10-7

6097.56

3.6×10-7

30

-7

Qzmax

Table 4: Numerical average Sherwood number, dimensionless maximum velocities and their

sensitivities to the grid refinement for the three test cases.

Test case 1

Level (n) 3

2 3

Ne

Sh

Gd , n SenSh

Qxmax

SenQGdmax,n

Qzmax

SenQGdmax,n

22500

-3.77

2.6 × 10−7

20.90

6.2 × 10 −7

48.21

8.5 ×10−7

4

40000

-15.91

4.8 × 10 −7

88.86

2.8 × 10 −7

556.75

4.8 × 10−7

6

90000

-53.55

2.4×10-7

383.06

8.7×10-7

5695.41

8.1×10-7

792 793

31

x

z

Highlights

794 795 796 797

• We propose a new benchmark semi-analytical solution for density-driven flow problem

798

• The benchmark deals with a square porous cavity salted on one of its vertical wall

799

• We develop a semi-analytical solution using the Fourier-Galerkin method

800

• Semi-analytical solution is obtained for very small values of diffusion coefficient

801

• We show the worthiness of the new semi-analytical solution in codes benchmarking

802

32