Mass-Conservative Finite Volume Methods for the Richards equation - - PowerPoint PPT Presentation

Mass-Conservative Finite Volume Methods for the Richards’ equation

Lecturer: G. Manzini1

1IMATI - CNR, Pavia

Siena 2006

(Joint collaboration with Stefano Ferraris, Univ. of Turin, Italy)

Richards’ equation: mixed θ−ψ formulation

• The solution ψ(x, t) satisfies

∂θ(ψ) ∂t − div

• K(ψ)∇(ψ + z)
• = s,

in Ω ×

• Ω ⊂

• ψ, hydraulic pressure head;
• θ(ψ), volumetric water content;
• K(ψ), hydraulic conductivity tensor;
• z, vertical coordinate (positive upward)
• s, production/consumption source term;

Richards’ equation: head formulation

• The solution ψ(x, t) satisfies

η(ψ)∂ψ ∂t − div

• K(ψ)∇(ψ + z)
• = s,

in Ω ×

where

• η(ψ) = ∂θ(ψ)/∂ψ is the general storage term.
• Remark: the two formulations are theoretically equivalent,

but leads to different numerical discretizations.

To complete the model. . .

. . . we must assign:

• a proper set of boundary conditions:

Dirichlet BCs: ψ = gD,

• n

ΓD ×

Neumann BCs: −n · K(ψ)∇(ψ + z) = gN,

• n

ΓN ×

where ∂Ω = ΓD ∪ ΓN, and gD and gN are, respectively, Dirichlet and Neumann boundary data;

• an initial solution:

ψ(x, t = 0) = ψ0(x) for x ∈ Ω,

• the characteristic relations for θ(ψ), K(ψ), and η(ψ) depending
• n soil nature: Van Genuchten, e.g. HYDRUS-2D, 1999,

Brooks-Corey, Haverkamp (1977), and many other models. . .

Extension to multi-layered soils

Multi-layered soils are taken into consideration by assuming that Ω be partitioned into a set of soil zones (sub-surface layers) Ωl, l = 1, . . . , L such that Ω =

• l=1,...,L

Ωl, showing different constitutive relations: θ(ψ)|Ωl = θl(ψ), K(ψ)|Ωl = Kl(ψ), η(ψ)|Ωl = ηl(ψ).

The finite volume framework

Let Th =

• T
• be a suitable triangulation of Ω;

Let us integrate over any control volume T ∈ Th to obtain

• T

∂θ(ψ) ∂t dV −

• T

div

• K(ψ)∇ (ψ + z)
• dV =
• T

s dV; then, apply the Gauss-Green divergence theorem

• T

∂θ(ψ) ∂t dV −

• ∂T

n ·

• K(ψ)∇ (ψ + z)
• dS =
• T

s dV, introduce cell averages and split face contributions to flux balance d dt

• T

θ(ψ) dV −

• e∈∂T

|e|

• ∂T

n ·

• K(ψ)∇ (ψ + z)
• dS =
• T

s dV.

Finally, let us introduce the cell-average approximation of the pressure field, e.g. ψh ≡ {ψT}, and the numerical flux, Gij(ψh) = |eij| nij · Kij

• G⋄

ij(ψh) + ζ

• ,

where ζ = ∇z is the versor along Z, 1 Kij(ψh) = 1 2

• 1

K(ψj) + 1 K(ψi)

• ,

and G⋄

ij(ψh) is the face-based gradient on the diamond Dij that is

given by G⋄

ij(ψh) = G(n) ij (ψh)nij + ψβ − ψα

|eij| tij.

The normal component G(n)

ij (ψh) of the face-based gradient G⋄ ij(ψh) is

G(n)

ij (ψh) =

1 |Dij|

• k

w(n)

ij,k ψk + g(n) ij

where the weights

• w(n)

ij,k

• are obtained by the Least Squares

reconstruction algorithm. Rewriting

• G(ψ)ψ + g(ψ)
• i =
• eij∈∂Ti

|eij| nij · KijG⋄

ij(ψh),

• h(ψ)
• i =
• eij∈∂Ti

|eij| nij · Kijζ, we finally obtain the vector formulation dθ(ψ) dt

• i −
• G(ψ)ψ + g(ψ) + h(ψ)
• i = si(ψ).

Time discretization

• Crucial issue: we are advancing in time θ(ψ), but we need ψ

for flux balancing: dθ(ψ) dt

• i
• non-linear

function of ψ +

• G(ψ)ψ + g(ψ)
• non-linear

flux term + h(ψ) gravity term

• i = si(ψ).
• the simplest approach is to switch to a discrete form of the

head-based formulation by using dθi(ψi) dt ≈ η(ψi)dψi dt .

• Using the implicit Euler F.D. scheme for the time derivative give

us η(ψn+1

i

)ψn+1

i

− ψn

i

∆t −

• G(ψn+1)ψn+1 + g(ψn+1) + h(ψn+1)
• i = s(ψn+1).
• This non-linear algebraic problem can be solved by using a

fixed-point iterative scheme (Picard):

• diag
• η(ψn+1,m

i

)

• − ∆tG(ψn+1,m)
• ψn+1,m+1 =

η(ψn+1,m)ψn + ∆t

• si(ψn+1,m) + g(ψn+1,m) + h(ψn+1,m)

• The (so-called) delta-form is obtained by solving for δm+1,m
• diag
• η(ψn+1,m)
• − ∆tG(ψn+1,m)

ψn+1,m+1 − ψn+1,m

• =δm+1,m

= +∆t

• s(ψn+1,m) + g(ψn+1,m) + h(ψn+1,m)
• −diag
• η(ψn+1,m)
• ψn+1,m − ψn

+ ∆tG(ψn+1,m)ψn+1,m

• Iterate on m to get ψn+1,m+1 from ψn+1,m starting from

ψn+1,0 = ψn until

• δm+1,m

≤ TOL and, finally, set ψn+1,m+1 = ψn+1,m + δm+1,m This approach is easy, but shows very poor mass conservation!

What is a mass-conservative discretization?

Following M. Celia et al., W.R.R. 1990, let us define:   total additional mass in the domain   =   total mass in the domain at any time t > 0   −   total mass in the domain at initial time t = 0   ,

• total net flux

into the domain

• =

flux balance integrated from initial time t = 0 to the current time t

• MASS BALANCE RATIO

= total additional mass in the domain total net flux into the domain Clearly, it must be MASS BALANCE RATIO=1.

How does the head-based approach perform?

Let us consider, for example, the vertical column infiltration problems proposed by M. Celia et al. W.R.R. 1990: Case 1 Data: column range: [bottom] 0 ≤ ζ ≤ 40 [top] bottom pressure: −61.5 cm top pressure: −20.7 cm initial pressure: −61.5 cm for 0 ≤ ζ < 40 running time: from T = 0 to T = 360 s. Constitutive relationships (Haverkamp et al.): θ(ψ) = α(θs − θr) α + |ψ|β + θr K(ψ) = Ks A A + |ψ|γ (α, β, γ, θs, θr, Ks and A are known parameters.)

Column infiltration problem: Case 1

100 200 300 400 Time-step size (sec) 0.75 0.8 0.85 0.9 0.95 1 Mass Balance Ratio Head-based 2D FV scheme (M. & Ferraris 2004) Head-based 1D FD scheme (Celia et.al., 1990)

Mass balance ratio for different time steps ∆t

How does the head-based approach perform?

Case 2 Data: column range: [bottom] 0 ≤ ζ ≤ 100 [top] bottom pressure: −1000 cm top pressure: −120 cm initial pressure: −1000 cm for 0 ≤ ζ < 100 running time: from T = 0 to T = 1 day. Constitutive relationships (Van Genuchten et al.): θ(ψ) = θs − θr 1 + (ε |ψ|n)m + θr K(ψ) = Ks

• 1 − (ε |ψ|)n−1

1 + (ε |ψ|n) −m2

• 1 + (ε |ψ|n)

m

2

(m, n, ε, θs, θr, and Ks are known parameters.)

Column infiltration problem: Case 2

1000 2000 3000 4000 Time-step size (sec) 0.4 0.5 0.6 0.7 0.8 0.9 1 Mass Balance Ratio Head-based 2D FV scheme (M. & Ferraris, 2004) Head-based 1D FD scheme (Celia et.al. 1990)

Mass balance ratio for different time steps ∆t

How can we improve this behavior?

• Main Idea: M. Celia et al. observed that mass conservation may

be lost in the head-based formulation because of a poor approximation of the time derivative of θ(ψ): ∂θ(ψ) ∂t

• n+1

i

≈ 1 |Ti|

• Ti

∂θ(ψ) ∂t

• n+1

i

dV ≈ d dt θ(ψi(t))

• n+1

≈ η(ψn

i )ψn+1 i

− ψn

i

∆t + O(|δψ|)

• Better strategy: develop θ(ψ) to second-order in ψ-variation in

the cell-average integral θn+1,m+1

i

= 1 |Ti|

• Ti

θ(ψn+1,m+1

i

) dV

How can we improve this behavior?

A straightforward calculation gives θn+1,m+1

i

= 1 |Ti|

• Ti

θ(ψn+1,m+1

i

) dV = 1 |Ti|

• Ti

θ

• ψn+1,m

i

+ (ψn+1,m+1

i

− ψn+1,m

i

)

• = δm+1,m

i

• dV

= 1 |Ti|

• Ti
• θ(ψn+1,m

i

) +

• ∂θ

∂ψ

• n+1,m

i

• ψn+1,m+1

i

− ψn+1,m

i

• + O(|ψ|2)
• dV

≈ θn+1,m

i

+ η(ψn+1,m)δm+1,m

i

+ O(|δψ|2).

• Let us start again from

dθ dt

• i −
• G(ψ)ψ + g(ψ) + h(ψ)
• i = si(ψ),
• and discretize directly the time derivative of θ(ψ) to get

θn+1

i

− θn

i

∆t −

• G(ψn+1)ψn+1 + g(ψn+1) + h(ψn+1)
• i = si(ψn+1).
• Use Picard iterations for solving the non-linear problem

θn+1,m+1

i

− θn

i

∆t −

• G(ψn+1,m)ψn+1,m+1 + g(ψn+1,m) + h(ψn+1,m)
• i = si(ψn+1,m).

• Re-write Picard iterations in δ-form,

θn+1,m+1

i

− θn+1,m

i

∆t −

• G(ψn+1,m)
• ψn+1,m+1 − ψn+1,m

+ g(ψn+1,m) + h(ψn+1,m)

• i =

si(ψn+1,m) + G(ψn+1,m)ψn+1,m

• i − θn+1,m

i

− θn

i

∆t , (recall δm+1,m = ψn+1,m+1 − ψn+1,m).

• We will use the previously derived second-order development

θn+1,m+1

i

≈ θn+1,m

i

+ η(ψn+1,m)δm+1,m

i

+ O(|δψ|2).

• Finally, we obtain the vector expression:
• diag
• η(ψn+1,m)
• − ∆tG(ψn+1,m)
• δm+1,m =

∆t

• s(ψn+1,m) + g(ψn+1,m) + h(ψn+1,m)
• −
• θn+1,m − θn

+ ∆tG(ψn+1,m)ψn+1,m,

• to be compared with the head-based discretization:
• diag
• η(ψn+1,m)
• − ∆tG(ψn+1,m)
• δm+1,m =

∆t

• s(ψn+1,m) + g(ψn+1,m) + h(ψn+1,m)
• −diag
• η(ψn+1,m)
• ψn+1,m − ψn

+ ∆tG(ψn+1,m)ψn+1,m.

Evaluation of mass conservation, Case 1

100 200 300 400

Time-step size (sec)

0.75 0.8 0.85 0.9 0.95 1

Mass Balance Ratio

Mixed-based (Celia et.al. 1990) 2D FV Standard Head-based 2D FV Standard Head-based 1D FD

Evaluation of mass conservation, Case 2

1000 2000 3000 4000

Time-step size (sec)

0.4 0.5 0.6 0.7 0.8 0.9 1

Mass Balance Ratio

Mixed-based (Celia et.al. 1990) 2D FV Standard Head-based 2D FV Standard Head-based 1D FD

Evaluation of mass conservation, Case 1

10 20 30 40

Z (cm)

• 60
• 50
• 40
• 30
• 20

Pressure Head (cm)

t=120s t=240s t=360s t=480s t=600s

Evaluation of mass conservation, Case 1

10 20 30 40

Z (cm)

0.1 0.15 0.2 0.25

Water Content

t=120s t=240s t=360s t=480s t=600s

Evaluation of mass conservation, Case 2

20 40 60 80 100

Z (cm)

• 1000
• 800
• 600
• 400
• 200

Pressure Head (cm)

t=12h t=24h t=36h t=48h

Evaluation of mass conservation, Case 2

20 40 60 80 100

Z (cm)

0.1 0.15 0.2

Water Content

t=12h t=24h t=36h t=48h

Multi-layered Calculations

The performance of the Finite Volume method for predicting unsaturated flow in multi-layered calculation is measured in two different test cases:

• horizontally layered soil (three layers), e.g. Baca, Chung,

Mulla, IJNMF, 1997;

• curvilinearly layered soil (two layers), M. & Ferraris, AdWR

2004.

Horizontally layered soil

Soil Data: range: 0cm ≤ ζ ≤ 25.5cm [top] surface crust: 25cm ≤ ζ ≤ 25.5cm Ks = 0.0616cm h−1 tilled layer: 15cm ≤ ζ ≤ 25cm Ks = 1.396cm h−1 sub-soil: 0cm ≤ ζ ≤ 15cm Ks = 0.312cm h−1 Constitutive relationships (Brooks-Corey et al.): θ(ψ) =

• θs

ψ ≥ ψa, θs(ψ/ψa)−1/b ψ < ψa K(ψ) =

• Ks

ψ ≥ ψa, Ks(ψ/ψa)−1/b ψ < ψa (θs, ψa, and b are known parameters.)

Simulation Data: Case 1 bottom pressure: −80 cm top pressure: 0 cm initial pressure: −80 cm for 0 ≤ ζ < 25.5 cm running time: from T = 0 to T = 4 hours. Case 2 bottom pressure: −1000 cm top pressure: 0 cm initial pressure: −1000 cm for 0 ≤ ζ < 25.5 cm running time: from T = 0 to T = 4 hours.

Horizontally layered soil, Case 1

5 10 15 20 25

Z (cm)

• 100
• 80
• 60
• 40
• 20

Pressure Head (cm)

t=30m t=1h t=1h30m t=2h

Horizontally layered soil, Case 1

5 10 15 20 25

Z (cm)

0.35 0.4 0.45 0.5 0.55

Water Content

t=0 t=1h t=2h t=3h t=4h t=1h t=2h t=3h t=4h

Horizontally layered soil, Case 2

5 10 15 20 25

Z (cm)

• 1000
• 800
• 600
• 400
• 200

Pressure Head (cm)

t=1h t=2h t=3h t=4h

Horizontally layered soil, Case 2

5 10 15 20 25

Z (cm)

• 50
• 40
• 30
• 20
• 10

Pressure Head (cm)

t=1h t=2h t=3h t=4h

Horizontally layered soil, Case 2

5 10 15 20 25

Z (cm)

0.2 0.3 0.4 0.5 0.6

Water Content

t=1h t=2h t=3h t=4h t=0 t=1h t=2h t=3h t=4h

Curvilinearly layered soil

Simulation Data domain: [0, Lx] × [0, Lz], with Lx = Lz100 cm bottom pressure: 0 cm top pressure: 0 cm initial pressure: (100 − z) cm for 0 ≤ z < 25.5 cm vertical BCs: homogeneous Neumann (no flux) running time: from T = 0 to T = 2 days. The curvilinear layer is the curve (x, ζ(x)) such that ζ(x) = 1 10

• 1 − cos(πx/Lx)
• + 0.45
• Lζ

Van Genuchten constitutive curves are used, in particular with layer 1: z ≥ ζ(x) Ks = 0.25, layer 2: z ≤ ζ(x) Ks = 2.00.

Curvilinearly layered soil

100 100

Curvilinearly layered soil

100 100 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 0.42 0.44 0.46 0.48

T = 0

100

100

0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 0.42 0.44 0.46 0.48

T = 0.25 day

Curvilinearly layered soil

100

100

0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 0.42 0.44 0.46 0.48

T = 0.5 day

100

100

0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 0.42 0.44 0.46 0.48

T = 0.75 day

Curvilinearly layered soil

100

100

0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 0.42 0.44 0.46 0.48

T = 1 day

100

100

0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 0.42 0.44 0.46 0.48

T = 1.25 day

Curvilinearly layered soil

100

100

0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 0.42 0.44 0.46 0.48

T = 1.5 day

100

100

0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 0.42 0.44 0.46 0.48

T = 1.75 day

Curvilinearly layered soil

100

100

0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 0.42 0.44 0.46 0.48

T = 2 day

Final Remarks

• The second-order Finite Volume method based on

diamond discretization of fluxes is valid alternative to more “traditional” numerical methods like FDM, FEM, for solving flow models in partially saturated soils;

• the mixed θ−ψ formulation of the Richards equation and

the mass-conservative technique by M. Celia et al. W.W.R 1990 makes it possible the development of a very accurate time-marching method;

• the scheme is naturally adapted to multi-layered soil

calculations.