Coupled Surface and Saturated/Unsaturated Ground Water Flow in - - PowerPoint PPT Presentation

Coupled Surface and Saturated/Unsaturated Ground Water Flow in Heterogeneous Media

Heiko Berninger∗, Ralf Kornhuber, and Oliver Sander Universit´ e de Gen` eve∗, Freie Universit¨ at Berlin and Matheon

Multiscale Simulation & Analysis in Energy and the Environment, Radon Special Semester 2011

Outline

Richards equation with homogeneous equations of state: Berninger, Kh. & Sander 10

• solver friendly finite element discretization: linear efficiency and robustness

Heterogeneous equations of state: Thesis of Berninger 07, Berninger, Kh. & Sander 07,09,...

• nonlinear domain decomposition: linear efficiency and robustness

Coupled Richards and shallow water equations: Ern et al. 06, Sochala et al. 09, Dawson 08, Berninger et al. 11

• continuous of mass flow and discontinuous pressure (clogging)
• mass conserving discretization (discontinuous Galerkin, ...) Dedner et al. 09, ...
• Steklov–Poincar´

e formulation and substructuring

• numerical experiments

All computations made with Dune Bastian, Gr¨

aser, Sander, ...

Outline

Richards equation with homogeneous equations of state: Berninger, Kh. & Sander 10

• solver friendly finite element discretization: linear efficiency and robustness

Heterogeneous equations of state: Thesis of Berninger 07, Berninger, Kh. & Sander 07,09,...

• nonlinear domain decomposition: linear efficiency and robustness

Coupled Richards and shallow water equations: Ern et al. 06, Sochala et al. 09, Dawson 08, Berninger et al. 11

• continuous of mass flow and discontinuous pressure (clogging)
• mass conserving discretization (discontinuous Galerkin, ...) Dedner et al. 09
• Steklov–Poincar´

e formulation and substructuring

• numerical experiments

All computations made with Dune Bastian, Gr¨

aser, Sander, ...

Runoff Generation for Lowland Areas

γE γSP saturated γD h γSP γE vadose

mathematical challenges:

• saturated/unsaturated ground water flow: non-smooth degenerate pdes

l Signorini-type bc (seepage face)

• coupling subsurface and surface water:

heterogeneous domain decomposition

• uncertain parameters (permeability, ...):

stochastic pdes Forster & Kh. 10, Forster 11

Saturated/Unsaturated Groundwater Flow: Richards Equation

∂ ∂t θ(p) + div v(x, p) = 0 , v(x, p) = −K(x) kr(θ(p))∇(p − ̺gz) equations of state: (Brooks & Corey, Burdine) θ(p) =

• θm + (θM − θm)
• p

pb

−ε (p ≤ pb) θM (p ≥ pb) kr(θ) = θ − θm θM − θm 3+2

ε

−5 −4 −3 −2 −1 1 2 3 4 5 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

saturation vs. pressure: p → θ(p)

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

relative permeability vs. saturation: θ → kr(θ)

Saturated/Unsaturated Groundwater Flow: Richards Equation

∂ ∂t θ(p) + div v(x, p) = 0 , v(x, p) = −K(x) kr(θ(p))∇(p − ̺gz) equations of state: (Brooks & Corey, Burdine) θ(p) =

• θm + (θM − θm)
• p

pb

−ε (p ≤ pb) θM (p ≥ pb) kr(θ) = θ − θm θM − θm 3+2

ε

−5 −4 −3 −2 −1 1 2 3 4 5 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

saturation vs. pressure: p → θ(p)

pb, ε → 0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

relative permeability vs. saturation: θ → kr(θ)

Homogeneous Equations of State Alt & Luckhaus 83, Otto 97

Kirchhoff Transformation: κ(p) := p kr(θ(q)) dq = ⇒ ∇κ(p) = kr(θ(p))∇p

−3 −1 3 −4/3 −1 3

generalized pressure: u := κ(p)

−2 −4/3 −1 2 0.21 0.95 1

M(u) := θ(κ−1(u)) separation of ill–conditioning and numerical solution: semilinear variational equation u(t) ∈ H1

0(Ω) :

• Ω

M(u)t v dx+

• Ω
• K∇u−kr(M(u))̺gez
• ∇v dx = 0

∀v ∈ H1

0(Ω)

Solver-Friendly Discretization

lumped implicit/explicit-upwind discretization in time, finite elements Sj ⊂ H1

0(Ω):

un+1

j

∈ Sj :

• Ω

ISj(M(un+1

j

) v) dx +

• Ω

τK∇un+1

j

∇v dx = ℓun

j (v)

∀v ∈ Sj equivalent convex minimization problem uj ∈ Sj : J (uj) + φj(uj) ≤ J (v) + φj(v) ∀v ∈ Sj quadratic energy J (v) = 1

2 (τK∇v, ∇v) − ℓun

j (v)

convex, l.s.c., proper functional φj(v) = Φ(v(p)) hp =

• Ω ISj(Φ(v)) dx

nonlinear convex function Φ : R → R ∪ {+∞} with ∂Φ = M

Algebraic Solution: Monotone Multigrid Kh. 99, 02, Gr¨

aser & Kh. 07

• given iterate uν

j

• fine grid smoothing:

− successive 1D minimization of J + φj in direction of nodal basis functions of Sj: 1 step of nonlinear Gauss–Seidel iteration → smoothed iterate ¯ uν

j

• coarse grid correction:

− Newton linearization of M(u) at ¯ uν

j

− constrain corrections to smooth regime of M 1 step of damped MMG → new iterate uν+1

j

= ⇒ (J + φj)(uν+1

j

) ≤ (J + φj)(uν

j)

0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8

• ¯

uν j (p) u M(u)

Theorem: (for V –cycle) Global convergence and asymptotic multigrid convergence ∀ −pb, ε ≥ 0 (robustness).

Algebraic Solution: Monotone Multigrid Kh. 99, 02, Gr¨

aser & Kh. 07

• given iterate uν

j

• fine grid smoothing:

− successive 1D minimization of J + φj in direction of nodal basis functions of Sj: 1 step of nonlinear Gauss–Seidel iteration → smoothed iterate ¯ uν

j

• coarse grid correction:

− Newton linearization of M(u) at ¯ uν

j

− constrain corrections to smooth regime of M 1 step of damped MMG → new iterate uν+1

j

= ⇒ (J + φj)(uν+1

j

) ≤ (J + φj)(uν

j)

0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8

• ¯

uν j (p) u M(u) Signorini condition

Theorem: (for V –cycle) Global convergence and asymptotic multigrid convergence ∀ −pb, ε ≥ 0 (robustness).

Solver–Friendly Discretization

• Kirchhoff transformation
• discretization: discrete minimization problem for uj
• algebraic solution by monotone multigrid (descent method!)
• inverse discrete Kirchhoff transformation

Solver–Friendly Discretization

• Kirchhoff transformation:

u = κ(p)

• discretization: discrete minimization problem for uj
• algebraic solution by monotone multigrid (descent method!)
• inverse discrete Kirchhoff transformation

Solver–Friendly Discretization

• Kirchhoff transformation:

u = κ(p)

• discretization: discrete minimization problem for uj
• algebraic solution by monotone multigrid (descent method!)
• inverse discrete Kirchhoff transformation

Solver–Friendly Discretization

• Kirchhoff transformation:

u = κ(p)

• discretization: discrete minimization problem for uj
• algebraic solution by monotone multigrid
• inverse discrete Kirchhoff transformation

Solver–Friendly Discretization

• Kirchhoff transformation:

u = κ(p)

• discretization: discrete minimization problem for uj
• algebraic solution by monotone multigrid
• discrete inverse Kirchhoff transformation:

pj = Ij(κ−1(uj))

Reinterpretation and Convergence Analysis Berninger, Kh. & Sander 10

reinterpretation in terms of physical variables: inexact finite element discretization with special quadrature points convergence properties: generalized variables: uj → u and M(uj) → M(u) in H1(Ω) physical variables pj → p and θj(pj) = Ij (M(uj)) → θ(p) in L2(Ω)

Reinterpretation and Convergence Analysis Berninger, Kh. & Sander 10

reinterpretation in terms of physical variables: inexact finite element discretization with special quadrature points convergence properties: generalized variables: uj → u and M(uj) → M(u) in H1(Ω) physical variables pj → p and θj(pj) = Ij (M(uj)) → θ(p) in L2(Ω)

Experimental Order of L2-Convergence

model problem: time discretized Richards equation without gravity physical parameters: Ω = (0, 2) × (0, 1), sandy soil → ε, θm, θM, pb, n triangulation; uniformly refined triangulation T11 (8 394 753 nodes)

1e-07 1e-06 1e-05 0.0001 0.001 0.01 0.1 1 0.0001 0.001 0.01 0.1 1 L2-error mesh size h 1e-07 1e-06 1e-05 0.0001 0.001 0.01 0.1 1 10 0.0001 0.001 0.01 0.1 1 L2-error mesh size h

generalized pressure u physical pressure p

Experimental Order of H1-Convergence

model problem: time discretized Richards equation without gravity physical parameters: Ω = (0, 2) × (0, 1), sandy soil → ε, θm, θM, pb, n triangulation; uniformly refined triangulation T11 (8 394 753 nodes)

0.001 0.01 0.1 1 10 0.0001 0.001 0.01 0.1 1 H1-error mesh size h 0.001 0.01 0.1 1 10 0.0001 0.001 0.01 0.1 1 H1-error mesh size h

generalized pressure u physical pressure p

Evolution of a Wetting Front in a Porous Dam

physical parameters: Ω = (0, 2) × (0, 1), sand → ε, θm, θM, pb, n triangulation; uniformly refined triangulation T4 (216 849 nodes) initial wetting front wetting front for t = 100s pressure pj

Efficiency and Robustness of Monotone Multigrid (MMG)

pre- and postsmoothing steps: V(3,3) cycle

20 40 60 80 100 0.2 0.4 0.6 0.8 1

from unsaturated... ...to saturated soil t = 0 convergence rates ρ over time t t = 250s

Robustness with respect to Soil Parameters

10

−2

10

−1

10 10

1

0.2 0.4 0.6 0.8 1 10

−2

10

−1

10 10

1

0.2 0.4 0.6 0.8 1

variation of ε: variation of −pb: ρave over ε ρave over −pb

Part II: Heterogeneous Equations of State

Spatially varying equations of state: no global Kirchhoff transformation! Ω = N

i=1 Ωi

θ = θ(x, p) = θi(p) ∀x ∈ Ωi kr = kr(x, θ) = kri(θ) ∀x ∈ Ωi

• Jumping soil parameters ε, pb

nonlinear domain decomposition:

Berninger & Kh. 07, Berninger 07, ...

physical transmission conditions: continuity of p and normal flux v · n generalized transmission conditions: κ−1

i (ui) = κ−1 j (uj) and ∇ui = ∇uj

nonlinear Dirichlet-Neumann and Robin methods efficient and robust subdomain solvers: monotone multigrid

Part II: Heterogeneous Equations of State

Spatially varying equations of state: no global Kirchhoff transformation! Ω = N

i=1 Ωi

θ = θ(x, p) = θi(p) ∀x ∈ Ωi kr = kr(x, θ) = kri(θ) ∀x ∈ Ωi

• Jumping soil parameters ε, pb

nonlinear domain decomposition:

Berninger & Kh. 07, Berninger 07, ...

physical transmission conditions: continuity of p and normal flux v · n generalized transmission conditions: κ−1

i (ui) = κ−1 j (uj) and ∇ui = ∇uj

nonlinear Dirichlet-Neumann and Robin methods efficient and robust subdomain solvers: monotone multigrid

Solver–Friendly (?) Discretization

• transmission problem
• local Kirchhoff transformation
• discretization: discrete nonlinear transmission problem
• algebraic solution by domain decomposition (fast local solvers!)
• local inverse discrete Kirchhoff transformation

Solver–Friendly (?) Discretization

• transmission problem
• local Kirchhoff transformation
• discretization: discrete nonlinear transmission problem
• algebraic solution by domain decomposition (fast local solvers!)
• local inverse discrete Kirchhoff transformation

Solver–Friendly (?) Discretization

• transmission problem
• local Kirchhoff transformation
• discretization: discrete nonlinear transmission problem
• algebraic solution by domain decomposition (fast local solvers!)
• local inverse discrete Kirchhoff transformation

Solver–Friendly (?) Discretization

• transmission problem
• local Kirchhoff transformation
• discretization: discrete nonlinear transmission problem
• algebraic solution by domain decomposition (fast local solvers!)
• local inverse discrete Kirchhoff transformation

Solver–Friendly (?) Discretization

• transmission problem
• local Kirchhoff transformation
• discretization: discrete nonlinear transmission problem
• algebraic solution by domain decomposition (fast local solvers!)
• local inverse discrete Kirchhoff transformation

Solver–Friendly (?) Discretization

• transmission problem
• local Kirchhoff transformation
• discretization: discrete nonlinear transmission problem
• algebraic solution by domain decomposition (fast local solvers!)
• local inverse discrete Kirchhoff transformation

Nonoverlapping Domain Decomposition Problem

Richards equation (time-discretized) in two different soils: θ(x, p) = θi(p), k(x, p) = ki(p) ∀x ∈ Ωi, i = 1, 2, θi monotonically increasing and Lipschitz, ki ∈ L∞(R), ki ≥ c > 0. Find p on Ω, p|Ωi = pi ∈ H1(Ωi), i = 1, 2, p|∂Ω = 0, such that θi(pi) − div(ki(pi)∇pi) = f

• n Ωi

p1 = p2

• n Γ

k1(p1)∇p1 · n = k2(p2)∇p2 · n

• n Γ

Ω1 Ω2 Γ n

Nonoverlapping Domain Decomposition Problem

Richards equation (time-discretized) in two different soils: θ(x, p) = θi(p), k(x, p) = ki(p) ∀x ∈ Ωi, i = 1, 2, θi monotonically increasing and Lipschitz, ki ∈ L∞(R), ki ≥ c > 0. Find p on Ω, p|Ωi = pi ∈ H1(Ωi), i = 1, 2, p|∂Ω = 0, such that θi(pi) − div(ki(pi)∇pi) = f

• n Ωi

p1 = p2

• n Γ

k1(p1)∇p1 · n = k2(p2)∇p2 · n

• n Γ

Ω1 Ω2 Γ n Kirchhoff transformation (different superposition operators on subdomains): κi(pi) = pi ki(q) dq = ui

Transformations in each Subdomain

Properties: − κ′

i = ki a.e. =

⇒ κi strictly increasing, Lipschitz and inverse Lipschitz − c2

Ωi |∇pi|2 ≤

• Ωi |∇ui|2 ≤ ki2

∞

• Ωi |∇pi|2 ⇒ ui|ΓH1/2

00 (Γ) ∼ pi|ΓH1/2 00 (Γ)

Transformed problem: Find u on Ω, u|Ωi = ui ∈ H1(Ωi), i = 1, 2, u|∂Ω = 0, such that Mi(ui) − ∆ui = f

• n Ωi

κ−1

1 (u1)

= κ−1

2 (u2)

• n Γ

∂u1 ∂n = ∂u2 ∂n

• n Γ

Analysis by Steklov–Poincar´ e Operators

• For given Dirichlet-value η ∈ Λ = H1/2

00 (Γ) find wi ∈ H1(Ωi), wi|∂Ω = 0, with

Mi(wi) − ∆wi = f in Ωi, wi = κi(η)

• n Γ

and set Siη = ∂wi

∂ni ∈ Λ′, the Neumann-value.

Analysis by Steklov–Poincar´ e Operators

• For given Dirichlet-value η ∈ Λ = H1/2

00 (Γ) find wi ∈ H1(Ωi), wi|∂Ω = 0, with

Mi(wi) − ∆wi = f in Ωi, wi = κi(η)

• n Γ

and set Siη = ∂wi

∂ni ∈ Λ′, the Neumann-value.

• With the Dirichlet-to-Neumann-maps or Steklov–Poincar´

e operators S1, S2: Find λ ∈ Λ such that Sλ = S1λ + S2λ = ∂u1 ∂n1 + ∂u2 ∂n2 = 0 in Λ′ is equivalent to the substructuring problem.

Nonlinear Dirichlet–Neumann Method

With damping parameter ϑ ∈ (0, 1) and given iterate λ0 on Γ :          M1(uk+1

1

) − ∆uk+1

1

= f

• n Ω1

uk+1

1

=

• n ∂Ω1 ∩ ∂Ω

uk+1

1

= κ1(λk)

• n Γ

         M2(uk+1

2

) − ∆uk+1

2

= f

• n Ω2

uk+1

2

=

• n ∂Ω2 ∩ ∂Ω

∂uk+1

2

∂n

=

∂uk+1

1

∂n

• n Γ

λk+1 := κ−1

2 (ϑ uk+1 2

|Γ + (1 − ϑ)κ2(λk))

Nonlinear Dirichlet–Neumann Method

With damping parameter ϑ ∈ (0, 1) and given iterate λ0 on Γ :          M1(uk+1

1

) − ∆uk+1

1

= f

• n Ω1

uk+1

1

=

• n ∂Ω1 ∩ ∂Ω

uk+1

1

= κ1(λk)

• n Γ

         M2(uk+1

2

) − ∆uk+1

2

= f

• n Ω2

uk+1

2

=

• n ∂Ω2 ∩ ∂Ω

∂uk+1

2

∂n

=

∂uk+1

1

∂n

• n Γ

λk+1 := κ−1

2 (ϑ uk+1 2

|Γ + (1 − ϑ)κ2(λk)) For Mi = 0 (stationary case): Convergence in 1D for small ϑ if 0 < c ≤ κ′

i(·) ≤ C.

Stationary Richards Equation: Strongly Heterogeneous Case

− domain: unit Yin Yang Ω = Ω1 ∪ Ω2 − hydrological data: USDA soil texture triangle (Maidment) Ω1 (grey clay): ε2 = 0.165, pb,2 = −0.373 [m], K1 = 1.67 · 10−7 [m/s] Ω2 (white sand): ε1 = 0.694, pb,1 = −0.073 [m], K2 = 6.54 · 10−5 [m/s] − Right hand side: source (f1 = 5 · 10−5) in grey circle sink (f2 = −2.5 · 10−3) in white circle − no gravity term Ω1 Ω2 − discretization: uniformly refined triangulations T0 – T7 (938 000 nodes)

Solution

Physical pressure p: Ranges: p1 ∈ [−56.1, 0.0] (clay) p2 ∈ [−36.2, 3.0] (sand) pi smooth across pi = pb,i p nonsmooth across Γ

Comparison of Nonlinear and Linear Dirichlet–Neumann Method

• linear case: mesh-independence for sufficiently small damping parameter (proof)
• nonlinear case:

(i) mesh-independence for sufficiently small damping parameter (proof in 1D) (ii) numerical mesh-independence (observed in 2D) (iii) mildly heterogeneous data: much damping (ϑopt ≈ 0.17) needed for acceptable optimal rates ρopt ≈ 0.77 (iv) strongly heterogeneous data: little damping (ϑopt ≈ 0.85) needed for good optimal rates ρopt ≈ 0.15

• general observation:

jumping diffusion coefficients K1/K2 ≫ 1 seem to improve convergence

Nonlinear Robin Method Gustave Robin (1855–1897)

With acceleration parameters γ1, γ2 > 0 and given u0

2 :

         M1(uk+1

1

) − ∆uk+1

1

= f

• n Ω1

uk+1

1

=

• n ∂Ω1 ∩ ∂Ω

∂uk+1

1

∂n

+ γ1κ−1

1 (uk+1 1

) =

∂uk

2

∂n + γ1κ−1 2 (uk 2)

• n Γ

         M2(uk+1

2

) − ∆uk+1

2

= f

• n Ω2

uk+1

2

=

• n ∂Ω2 ∩ ∂Ω

∂uk+1

2

∂n

− γ2κ−1

2 (uk+1 2

) =

∂uk+1

1

∂n

− γ2κ−1

1 (uk+1 1

)

• n Γ

Nonlinear Robin Method Gustave Robin (1855–1897)

With acceleration parameters γ1, γ2 > 0 and given u0

2 :

         M1(uk+1

1

) − ∆uk+1

1

= f

• n Ω1

uk+1

1

=

• n ∂Ω1 ∩ ∂Ω

∂uk+1

1

∂n

+ γ1κ−1

1 (uk+1 1

) =

∂uk

2

∂n + γ1κ−1 2 (uk 2)

• n Γ

         M2(uk+1

2

) − ∆uk+1

2

= f

• n Ω2

uk+1

2

=

• n ∂Ω2 ∩ ∂Ω

∂uk+1

2

∂n

− γ2κ−1

2 (uk+1 2

) =

∂uk+1

1

∂n

− γ2κ−1

1 (uk+1 1

)

• n Γ

convergence in 1D if 0 < c ≤ κ′

i(·) ≤ C and

Mi : R → R Lipschitz continuous and monotonically increasing, i = 1, 2.

Properties of the Linear Robin Method Lions 90, Gander 04, Liu 06, Dubois 07

• no mesh-independence in general. O.K.
• convergence speed can be increased by a good choice of Robin parameters. O.K.
• for γ1 = γ2 = γ one has the asymptotic behaviour

γopt = O(h−1/2) and ρopt = 1 − O(h1/2). O.K.

• for possibly different γ1 and γ2 one has the improved asymptotic behaviour

γopt = O(h−1/4) and ρopt = 1 − O(h1/4). quite O.K.

• for possibly different γ1 and γ2 and jumping diffusion coefficients

K1/K2 > 1

• ne can even obtain

ρopt = 1/µ − O(h1/4/µ), i.e. mesh-independence. O.K.

• optimal Robin tends to undamped Dirichlet-Neumann for

K1/K2 → ∞ O.K.

Comparison of Nonlinear and Linear Robin Method

• no mesh-independence in general. O.K.
• convergence speed can be increased by a good choice of Robin parameters. O.K.
• for γ1 = γ2 = γ one has the asymptotic behaviour

γopt = O(h−1/2) and ρopt = 1 − O(h1/2). O.K.

• for possibly different γ1 and γ2 one has the improved asymptotic behaviour

γopt = O(h−1/4) and ρopt = 1 − O(h1/4). quite O.K.

• for possibly different γ1 and γ2 and jumping diffusion coefficients

K1/K2 > 1

• ne can even obtain

ρopt = 1/µ − O(h1/4/µ), i.e. mesh-independence. O.K.

• optimal Robin tends to undamped Dirichlet-Neumann for

K1/K2 → ∞ O.K.

Robustness: Dirichlet–Neumann Berninger, Kh. & Sander 11

problem: unit Yin Yang problem with homogeneous soil (sand / sand)

10

−2

10

−1

10 10

1

0.2 0.4 0.6 0.8 1 −10

1

−10 −10

−1

−10

−2

0.2 0.4 0.6 0.8 1

variation of ε: variation of pb: ρave over ε ρave over pb

Conclusion of Numerical Experiments Berninger, Kh. & Sander 09 Nonlinear Dirichlet-Neumann and Robin methods behave surprisingly similar as in the linear case.

Part III: Coupling of Ground and Surface Water

non-moving surface water: compartment model mass conservation: m′(t) =

• γD∪γSP v(x, p) · n dσ

hydrostatic pressure: p = h(m)ρg on γD, Signorini b.c. on γSP ∪ γE

γE γSP saturated γD h γSP γE vadose

discretization: explicit Euler method

Compartment Model with Heterogeneous Dam Berninger 07

flow of water through into a 2D-domain with 4 different soils and surface water: − soil types (from top to bottom): sand, loamy sand, sandy loam, loam soil parameters: εi ∈ [0.252, 0.694], Ki ∈ [3.67 · 10−6, 6.54 · 10−5], pb,i ∈ [−0.147, −0.073], i = 1, 2, 3, 4 − initial condition: p0 = −10 (dry soil) and surface water − boundary conditions: dynamic Dirichlet/Signorini and surface water on the top −v · n = 3 · 10−4 [m/s] on the left and v · n = 0 elsewhere − discretization: implicit/explicit Euler (time step size 10 [s]), FE (h = 1/24) − domain decomposition: nonlinear Robin method (γ = 10−4) − local solver: monotone multigrid (3 levels, accuracy 10−12 in H1-norm w.r.t. u)

Coupling of Richards Equation with the Shallow Water Equations

(R) θ(p)t + div v(x, p) = 0 , v(x, p) = −K(x) kr(θ(p))∇(p − ̺gz)

• n Ω ⊂ Rd

(S1) mass conservation: ht + div q = gS

• n Γ

(S2) momentum conservation: qt + div (q2/h + gh2/2) = 0

• n Γ

coupling: mass flux gS = v · n hydrostatic pressure: p = ̺gh

• n Γ

Γ Ω

Coupling of Richards Equation with the Shallow Water Equations

(R) θ(p)t + div v(x, p) = 0 , v(x, p) = −K(x) kr(θ(p))∇(p − ̺gz)

• n Ω ⊂ Rd

(S1) mass conservation: ht + div q = gS

• n Γ

(S2) momentum conservation: qt + div (q2/h + gh2/2) = 0

• n Γ

coupling: mass flux: gS = v · n hydrostatic pressure: p = ̺gh

• n Γ

Γ Ω

Coupling of Richards Equation with the Shallow Water Equations

(R) θ(p)t + div v(x, p) = 0 , v(x, p) = −K(x) kr(θ(p))∇(p − ̺gz)

• n Ω ⊂ Rd

(S1) mass conservation: ht + div q = gS

• n Γ

(S2) momentum conservation: qt + div (q2/h + gh2/2) = 0

• n Γ

coupling: mass flux: gS = v · n hydrostatic pressure: p = ̺gh

• n Γ

Γ Ω

Coupling of Richards Equation with the Shallow Water Equations

(R) θ(p)t + div v(x, p) = 0 , v(x, p) = −K(x) kr(θ(p))∇(p − ̺gz)

• n Ω ⊂ Rd

(S1) mass conservation: ht + div q = gS

• n Γ

(S2) momentum conservation: qt + div (q2/h + gh2/2) = 0

• n Γ

coupling: mass flux: gS = v·n

• discont. pressure:

p = ̺gh+α div(q/h) on Γ J¨ ager–Mikeli´ c condition (Navier–Stokes/Darcy): pD = pN − α n · ∂uN

∂n

• n Γ

and n · ∂uN

∂n = O(ε), ε = char. pore size.

shallow water equations on Γ: n · ∂uN

∂n ≈ ∂ ∂zuz = − ∂ ∂xux − ∂ ∂yuy = − div(q/h).

Coupling of Richards Equation with the Shallow Water Equations

(R) θ(p)t + div v(x, p) = 0 , v(x, p) = −K(x) kr(θ(p))∇(p − ̺gz)

• n Ω ⊂ Rd

(S1) mass conservation: ht + div q = gS

• n Γ

(S2) momentum conservation: qt + div (q2/h + gh2/2) = 0

• n Γ

coupling: mass flux: gS = v · n

• discont. pressure:

p = ̺gh + αv · n

• n Γ

clogging: pressure discontinuity due to a thin and nearly impermeable river bed with thickness ε and conductivity Kε. Darcy’s law provides: v = −Kε∇peff := −Kε ̺gh − p ε n , α = ε Kε (leakage coefficient)

Formal Steklov–Poincar´ e Formulation Quarteroni & Valli 99, Deparis et al. 05, ...

coupled problem:

Sochala et al. 09, Dawson 08, ...

n θ(p)t + div v(x, p) = 0 , v(x, p) = −K(x) kr(θ(p))∇(p − ̺gz) on Ω × [0, T] ht + div q = gS, qt + div (q2/h + gh2/2) = 0

• n Γ × [0, T]

gS = v · n, p = ̺gh

• n Γ × [0, T]

Dirichlet to Neumann maps: subsurface: h → SΩ(h) = v · n surface: h → SΓ(h) = − div q(h) Steklov–Poincar´ e formulation: ht = SΩ(h) + SΓ(h)

• n Γ × [0, T]

Formal Steklov–Poincar´ e Formulation Quarteroni & Valli 99, Deparis et al. 05, ...

coupled problem:

Sochala et al. 09, Dawson 08, ...

n θ(p)t + div v(x, p) = 0 , v(x, p) = −K(x) kr(θ(p))∇(p − ̺gz) on Ω × [0, T] ht + div q = gS, qt + div (q2/h + gh2/2) = 0

• n Γ × [0, T]

gS = v · n, p = ̺gh

• n Γ × [0, T]

Dirichlet to Neumann maps: subsurface: h → SΩ(h) = v · n surface: h → SΓ(h) = − div q(h) Steklov–Poincar´ e formulation: ht = SΩ(h) + SΓ(h)

• n Γ × [0, T]

Formal Steklov–Poincar´ e Formulation Quarteroni & Valli 99, Deparis et al. 05, ...

coupled problem:

Sochala et al. 09, Dawson 08, ...

n θ(p)t + div v(x, p) = 0 , v(x, p) = −K(x) kr(θ(p))∇(p − ̺gz) on Ω × [0, T] ht + div q = gS, qt + div (q2/h + gh2/2) = 0

• n Γ × [0, T]

gS = v · n, p = ̺gh

• n Γ × [0, T]

Dirichlet to Neumann maps: subsurface: h → SΩ(h) = v · n surface: h → SΓ(h) = − div q(h) Steklov–Poincar´ e formulation: ht = SΩ(h) + SΓ(h)

• n Γ × [0, T]

Formal Steklov–Poincar´ e Formulation Quarteroni & Valli 99, Deparis et al. 05, ...

coupled problem:

Sochala et al. 09, Dawson 08, ...

n θ(p)t + div v(x, p) = 0 , v(x, p) = −K(x) kr(θ(p))∇(p − ̺gz) on Ω × [0, T] ht + div q = gS, qt + div (q2/h + gh2/2) = 0

• n Γ × [0, T]

gS = v · n, p = ̺gh

• n Γ × [0, T]

Dirichlet to Neumann maps: subsurface: h → SΩ(h) = v · n surface: h → SΓ(h) = − div q(h) Steklov–Poincar´ e formulation: ht = SΩ(h) + SΓ(h)

• n Γ × [0, T]

Iterative Substructuring

Dirichlet–Neumann iteration: hν+1

t

− SΓ(hν+1) = SΩ(hν),

• n Γ × [0, T],

k = 1, . . . , N local Dirichlet–Neumann iteration: hν+1

k,t − SΓ(hν+1 k

) = SΩ(hν

k),

h0

k = hν∗ k−1

• n Γ × [tk−1, tk],

k = 1, . . . , N Dirichlet–Neumann time stepping: hk,t + SΓ(hk) = SΩ(hk−1),

• n Γ × [tk−1, tk],

k = 1, . . . , N

• ther substructuring methods:

Neumann–Dirichlet, Neumann–Neumann, ... Neumann–Dirichlet: hν+1

t

− SΩ(hν+1) = SΓ(hν), Neumann–Neumann, ... discretization: different schemes, time steps, grids, ...

Iterative Substructuring

Dirichlet–Neumann iteration: hν+1

t

− SΓ(hν+1) = SΩ(hν),

• n Γ × [0, T],

k = 1, . . . , N local Dirichlet–Neumann iteration: hν+1

k,t − SΓ(hν+1 k

) = SΩ(hν

k),

h0

k = hν∗ k−1

• n Γ × [tk−1, tk],

k = 1, . . . , N Dirichlet–Neumann time stepping: hk,t + SΓ(hk) = SΩ(hk−1),

• n Γ × [tk−1, tk],

k = 1, . . . , N

• ther substructuring methods:

Neumann–Dirichlet, Neumann–Neumann, ... Neumann–Dirichlet: hν+1

t

− SΩ(hν+1) = SΓ(hν), Neumann–Neumann, ... discretization: different schemes, time steps, grids, ...

Iterative Substructuring

Dirichlet–Neumann iteration: hν+1

t

− SΓ(hν+1) = SΩ(hν),

• n Γ × [0, T],

k = 1, . . . , N local Dirichlet–Neumann iteration: hν+1

k,t − SΓ(hν+1 k

) = SΩ(hν

k),

h0

k = hν∗ k−1

• n Γ × [tk−1, tk],

k = 1, . . . , N Dirichlet–Neumann time stepping: hk,t − SΓ(hk) = SΩ(hk−1),

• n Γ × [tk−1, tk],

k = 1, . . . , N

• ther substructuring methods:

Neumann–Dirichlet, Neumann–Neumann, ... Neumann–Dirichlet: hν+1

t

− SΩ(hν+1) = SΓ(hν), Neumann–Neumann, ... discretization: different schemes, time steps, grids, ...

Iterative Substructuring

Dirichlet–Neumann iteration: hν+1

t

− SΓ(hν+1) = SΩ(hν),

• n Γ × [0, T],

k = 1, . . . , N local Dirichlet–Neumann iteration: hν+1

k,t − SΓ(hν+1 k

) = SΩ(hν

k),

h0

k = hν∗ k−1

• n Γ × [tk−1, tk],

k = 1, . . . , N Dirichlet–Neumann time stepping: hk,t − SΓ(hk) = SΩ(hk−1),

• n Γ × [tk−1, tk],

k = 1, . . . , N

• ther substructuring methods:

Neumann–Dirichlet: hν+1 = S−1

Ω (hν t − SΓ(hν)),

Neumann–Neumann, ... discretization: different schemes, time steps, grids, ...

Iterative Substructuring

Dirichlet–Neumann iteration: hν+1

t

− SΓ(hν+1) = SΩ(hν),

• n Γ × [0, T],

k = 1, . . . , N local Dirichlet–Neumann iteration: hν+1

k,t − SΓ(hν+1 k

) = SΩ(hν

k),

h0

k = hν∗ k−1

• n Γ × [tk−1, tk],

k = 1, . . . , N Dirichlet–Neumann time stepping: hk,t − SΓ(hk) = SΩ(hk−1),

• n Γ × [tk−1, tk],

k = 1, . . . , N

• ther substructuring methods:

Neumann–Dirichlet: hν+1 = S−1

Ω (hν t − SΓ(hν)),

Neumann–Neumann, ... discretization: different schemes, time steps, grids, ...

Coupling Conditions Revisited

continuous pressure: Dirichlet conditions: pν+1

Γ

(t) = ρghν(t) regularity gap: hν(t) ∈ L∞(Γ) ⊂ H1/2(Γ) clogging: p = ρgh + αv · n Robin condition in weak form:

• Γ

(pν+1 − ρghν)v dσ well-defined: hν(t) ∈ L1(Γ) ∩ L∞(Γ) ⊂ H1/2(Γ)′ = ⇒ Robin-Neumann iteration

Coupling Conditions Revisited

continuous pressure: Dirichlet conditions: pν+1

Γ

(t) = ρghν(t) regularity gap: hν(t) ∈ L∞(Γ) ⊂ H1/2(Γ) = ⇒ pν+1(t) ∈ H1(Ω) clogging: p = ρgh + αv · n Robin condition in weak form:

• Γ

(pν+1 − ρghν)v dσ well-defined: hν(t) ∈ L1(Γ) ∩ L∞(Γ) ⊂ H1/2(Γ)′ = ⇒ Robin-Neumann iteration

Coupling Conditions Revisited

continuous pressure: Dirichlet conditions: pν+1

Γ

(t) = ρghν(t) regularity gap: hν(t) ∈ L∞(Γ) ⊂ H1/2(Γ) = ⇒ pν+1(t) ∈ H1(Ω) discontinuous pressure (clogging): pν+1 = ρghν + αvν+1 · n Robin condition in weak form:

• Γ

(pν+1 − ρghν)v dσ well-defined: hν(t) ∈ L∞(Γ) ⊂ H1/2(Γ)′ = ⇒ Robin-Neumann iteration

Coupling Conditions Revisited

continuous pressure: Dirichlet conditions: pν+1

Γ

(t) = ρghν(t) regularity gap: hν(t) ∈ L∞(Γ) ⊂ H1/2(Γ) = ⇒ pν+1(t) ∈ H1(Ω) discontinuous pressure (clogging): pν+1 = ρghν + αvν+1 · n weak formulation of Robin condition: α−1 pν+1, vΓ − (ρghν, v)

• Γ

well-defined: hν(t) ∈ L∞(Γ) ⊂ H1/2(Γ)′ = ⇒ Robin-Neumann iteration

Coupling Conditions Revisited

continuous pressure: Dirichlet conditions: pν+1

Γ

(t) = ρghν(t) regularity gap: hν(t) ∈ L∞(Γ) ⊂ H1/2(Γ) = ⇒ pν+1(t) ∈ H1(Ω) discontinuous pressure (clogging): pν+1 = ρghν + αvν+1 · n weak formulation of Robin condition: α−1 pν+1, vΓ − (ρghν, v)

• Γ

well-defined: hν(t) ∈ L∞(Γ) ⊂ H1/2(Γ)′ = ⇒ Robin-Neumann iteration

Coupling Conditions Revisited

continuous pressure: Dirichlet conditions: pν+1

Γ

(t) = ρghν(t) regularity gap: hν(t) ∈ L∞(Γ) ⊂ H1/2(Γ) = ⇒ pν+1(t) ∈ H1(Ω) discontinuous pressure (clogging): pν+1 = ρghν + αvν+1 · n weak formulation of Robin condition: α−1 pν+1, vΓ − (ρghν, v)

• Γ

well-defined: hν(t) ∈ L∞(Γ) ⊂ H1/2(Γ)′ Robin-Neumann iteration: hν+1

t

− SΓ(hν+1) = SΩ(hν),

Discretization of the Shallow Water Equations

spatial discretization of :

∂ ∂tu = S − div F(u):

uh ∈ VΓ :

∂ ∂t

• Γ

uhw dσ =

N

• i=1
• Γi

S w − F(uh) · ∇w dσ −

• ∂Γi

w G(u+

h , u− h ) · n dσ

∀w ∈ V2

Γ

discontinuous Galerkin: VΓ = {v ∈ L2(Γ) | v|Γi polynomial of degree p ≥ 0} discrete flux G(u+

h , u− h ),

u+

h , u− h ∈ V2 Γ

time discretization by a stabilized Runge–Kutta scheme Dedner & Kl¨

• fkorn 08

Discrete Robin–Neumann Iteration

Robin boundary conditions for Richards equation: evaluation of right hand side:

• Γ

hν v dσ, hν ∈ VΓ, v ∈ Sh source term for shallow water equation: extension by zero: E : Sh|′

Γ → V′ Γ

E(v · n), w = v · n, wS, hierarchical decomposition w = wS + wDG discrete mass conservation:

• Ω

θν+1

h

dx +

• Γ

hν+1

h

dσ =

• Ω

θν

h dx +

• Γ

hν

h dσ ,

ν = 0, 1, . . . ,

Discrete Robin–Neumann Iteration

Robin conditions for Richards equation: ... + α−1pν+1, vΓ = α−1(hν, v)Γ, hν ∈ VΓ, v ∈ Sh source term for shallow water equation: (v · n, v)Γ = α−1(pν+1, v)Γ − (ρghν, v)Γ, pν+1

h

∈ Sh, v ∈ VΓ Discrete mass conservation:

• Ω

θν+1

h

dx +

• Γ

hν+1

h

dσ =

• Ω

θν

h dx +

• Γ

hν

h dσ ,

ν = 0, 1, . . . ,

Discrete Robin–Neumann Iteration

Robin conditions for Richards equation: ... + α−1pν+1, vΓ = α−1(hν, v)Γ, hν ∈ VΓ, v ∈ Sh source term for shallow water equation: (v · n, v)Γ = α−1 (pν+1, v)Γ − (ρghν, v)Γ

• ,

pν+1 ∈ Sh, v ∈ VΓ Discrete mass conservation:

• Ω

θν+1

h

dx +

• Γ

hν+1

h

dσ =

• Ω

θν

h dx +

• Γ

hν

h dσ ,

ν = 0, 1, . . . ,

Discrete Robin–Neumann Iteration

Robin conditions for Richards equation: ... + α−1pν+1, vΓ = α−1(hν, v)Γ, hν ∈ VΓ, v ∈ Sh source term for shallow water equation: (v · n, v)Γ = α−1 (pν+1, v)Γ − (ρghν, v)Γ

• ,

pν+1 ∈ Sh, v ∈ VΓ Discrete mass conservation:

• Ω

θν+1 dx +

• Γ

hν+1 dσ =

• Ω

θν dx +

• Γ

hν dσ + inflow − outflow

Numerical Results: Supercritical Surface Flow over Dry Soil

soil:

sandy soil, α = ρgL−1, L = 10−5s−1 (Wiese & N¨

utzmann 09)

problem:

p(0) = −20 Pa, h(0) = 10 m, q(0) = 10 m2s−1 + oscillating bc

discretization: ∆TΩ = 500s, ∆xΩ = 10·2−6m,

δTΓ = 3−1 ·10−5∆T , δxΓ = 10·400−1m

total mass conservation up to 10−10

Convergence of Robin–Neumann Iteration

α = ρg 105: weak coupling

Wiese & N¨ utzmann 09

α = ρg 5: medium coupling

Convergence of Robin–Neumann Iteration

α = ρg 105: weak coupling

Wiese & N¨ utzmann 09

α = ρg 5: medium coupling

Ongoing Work

Further convergence studies of Robin–Neumann iterations Other substructuring methods (Neumann–Robin, ...) Analysis of coupled problem (existence, uniqueness, ...) Numerical analysis of Robin–Neumann iterations (convergence, convergence rates, ...) flooding