Poster_density_driven_with_fracture_MLMC.pdf

Poster to be presented at Stochastic Numerics and Statistical Learning: Theory and Applications Workshop 2024, Kaust, Saudi Arabia, https://cemse.kaust.edu.sa/stochnum/events/event/snsl-workshop-2024.
In this work we have considered a setting that mimics the Henry problem \cite{Simpson2003,Simpson04_Henry}, modeling seawater intrusion into a 2D coastal aquifer. The pure water recharge from the ``land side'' resists the salinisation of the aquifer due to the influx of saline water through the ``sea side'', thereby achieving some equilibrium in the salt concentration. In our setting, following \cite{GRILLO2010}, we consider a fracture on the sea side that significantly increases the permeability of the porous medium.

The flow and transport essentially depend on the geological parameters of the porous medium, including the fracture. We investigated the effects of various uncertainties on saltwater intrusion. We assumed uncertainties in the fracture width, the porosity of the bulk medium, its permeability and the pure water recharge from the land side. The porosity and permeability were modeled by random fields, the recharge by a random but periodic intensity and the thickness by a random variable. We calculated the mean and variance of the salt mass fraction, which is also uncertain.

The main question we investigated in this work was how well the MLMC method can be used to compute statistics of different QoIs. We found that the answer depends on the choice of the QoI. First, not every QoI requires a hierarchy of meshes and MLMC. Second, MLMC requires stable convergence rates for $\EXP{g_{\ell} - g_{\ell-1}}$ and $\Var{g_{\ell} - g_{\ell-1}}$. These rates should be independent of $\ell$. If these convergence rates vary for different $\ell$, then it will be hard to estimate $L$ and $m_{\ell}$, and MLMC will either not work or be suboptimal. We were not able to get stable convergence rates for all levels $\ell=1,\ldots,5$ when the QoI was an integral as in \eqref{eq:integral_box}. We found that for $\ell=1,\ldots 4$ and $\ell=5$ the rate $\alpha$ was different. Further investigation is needed to find the reason for this. Another difficulty is the dependence on time, i.e. the number of levels $L$ and the number of sums $m_{\ell}$ depend on $t$. At the beginning the variability is small, then it increases, and after the process of mixing salt and fresh water has stopped, the variance decreases again.

The number of random samples required at each level was estimated by calculating the decay of the variances and the computational cost for each level. These estimates depend on the minimisation function in the MLMC algorithm.

To achieve the efficiency of the MLMC approach presented in this work, it is essential that the complexity of the numerical solution of each random realisation is proportional to the number of grid vertices on the grid levels.

1. Densitydrivenflowinfracturedporousmedia.UncertaintyquantificationviaMLMC.
AlexanderLitvinenko
1
,D.Logashenko
2
,R.Tempone
1,2
,G.Wittum
2
1
RWTHAachen,Germany,
2
KAUST,SaudiArabia
[email protected]
Abstract
Problem:Henry-likesaltwaterintrusioninfracturedmedia(nonlinear, time-
dependent,two-phaseflow)
Inputuncertainty:fracturewidth,porosity,permeability,andrecharge(mod-
elledbyrandomfields)
Solution:thesaltmassfraction(uncertainandtime-dependent)
Method:MultiLevelMonteCarlo(MLMC)method
Deterministicsolver:parallelmultigridsolverug4
Questions:
1.
2.
3.
variations?
4.
5.
time?ˆqin c= 0 c= 1 p=−ρ1gy 0 −1m 2m y x (2,−0.5) (1,−0.7) 1 2 3 4 5 6
(left)Henryproblemwithafracture.
(right)Massfractioncm(redforcm=1)
1.Henry-likeproblemsettings
Denote: the whole domainD, porous mediumM⊂ D, fracture porosityϕmand permeabilityKm, salt mass
fractioncm(t,x)andpressurepm(t,x). TheflowinMsatisfiesconservationlaws:
∂t(ϕmρm)+∇·(ρmqm)=0
∂t(ϕmρmcm)+∇·(ρmcmqm−ρmDm∇cm)=0
∑
x∈M (1)
withtheDarcy’slawforthevelocity:
qm=−
Km
µm
(∇pm−ρmg), x∈M, (2)
whereρm=ρ(cm)andµm=µ(cm)indicatethedensityandtheviscosityoftheliquidphase,Dm(t,x)denotesthe
moleculardiffusionandmechanicaldispersiontensor.
Thefractureisassumedtobefilledwithaporousmedium,too. Thefractureisrepresentedbyasurface S⊂D,
M∪S=D,M∩S=∅.
∂t(ϕ
fϵρ
f)+∇
S
·(ϵρ
fq
f)+Q
(1)
fn
+Q
(2)
fn
=0
∂t(ϕ
fϵρ
fc
f)+∇
S
·(ϵρ
fc
fq
f−ϵρ
fD
f∇
S
c
f)+P
(1)
fn
+P
(2)
fn
=0



x∈S, (3)
whereϵisthefracturewidth. TheDarcyvelocityalongthefractureis
q
f=−
K
f
µ
f
(∇
S
p
f−ρ
fg), x∈S. (4)
ThetermsQ
(k)
fn
andP
(k)
fn
,k∈{1,2},arethemassfluxesthroughthefacesS
(k)
ofthefracture:
Q
(k)
fn
:=ρ(c
(k)
m)q
(k)
fn
P
(k)
fn
:=ρ(c
(k)
m)c
(k)
upwind
q
(k)
fn
−ρ(c
(k)
m)D
(k)
fn
c
(k)
m−c
f
ϵ/2







x∈S
(k)
, (5)
Theuncertainwidthofthefracture,therecharge,andtheporosity ξ1,ξ2,ξ3∈U[−1,1],
ϵ(ξ1)=0.01·((1−0.01)·ξ1+(1+0.01))/2, (6)
ˆqx(t,ξ3)=3.3·10
−6
·(1+0.1ξ3)(1+0.1sin(πt/40)), (7)
ϕm(x,y,ξ2)=0.35·(1+0.02·(ξ2cos(πx/2)+ξ2sin(2πy)). (8)
Tocompute:cm.
Comput. domain:D×[0,T]. Wesetρ(c)=ρ0+(ρ1−ρ0)c,andD=ϕDI,
I.C.:c|t=0=0, c|x=2=1,p|x=2=−ρ1gy.c|x=0=0,ρq·ex|x=0=ˆqin.
Methods:Newton method, BiCGStab, preconditioned with the geometric multigrid method (V-cycle), ILU
β-
smoothersandGaussianelimination.
MultiLevelMonteCarlo(MLMC)method
SpatialandtemporalgridhierarchiesD0,D1,...,D
L,andT0,T1,...,T
L;n0=512,n
ℓ≈n0·4
ℓ
,τ
ℓ+1=
1
2
τ
ℓ,number
oftimestepsr
ℓ+1=2r
ℓorr
ℓ=r02
ℓ
. Computationcomplexityiss
ℓ=O(n
ℓr
ℓ)=O(2
3ℓγ
n0·r0)
MLMCapproximatesE[g
L]≈E[g]usingthefollowingtelescopicsum:
E[g
L]≈m
−1
0
m0X
i=1
g
(0,i)
0
+
L
X
ℓ=1

m
−1
ℓ
mℓX
i=1
(g
(ℓ,i)
ℓ
−g
(ℓ,i)
ℓ−1
)

.
MinimizeF(m0,...,m
L):=
P
L
ℓ=0
m
ℓs
ℓ+µ
2Vℓ
mℓ
,obtainm
ℓ=ε
−2
q
Vℓ
sℓ
P
L
i=0
√
Visi.
2.Numerics
Figure2:ThemeanvalueE[cm(t,x)]att={7τ,19τ,40τ,94τ}. Inallcases,E[cm](t,x)∈[0,1].
Figure3:ThevarianceV[cm(t,x)]att={7τ,19τ,40τ,94τ}. Maximalvalues(darkredcolour)ofV[cm]are1.9·10
−3
,
3.4·10
−3
,2.9·10
−3
,2.4·10
−3
respectively. Thedarkbluecolourcorrespondstoazerovalue.g 0
g 1
- g
0
g 2
- g
1
g 3
- g
2
g 4
- g
3
-6
-4
-2
0 g 0
g 1
- g
0
g 2
- g
1
g 3
- g
2
g 4
- g
3
-16
-14
-12
-10
-8
-6
-4 0 1 2 3 4
10
-3
10
-2 0 1 2 3 4
10
-8
10
-7
10
-6
10
-5
10
-4
Figure 4:(1-2)The weak and the strong convergences, QoI isg:=cm(t15,x1),α= 1.07,ζ1=−1.1,β= 1.97,
ζ2=−8. (3-4) Decay comparison of (left)E[g
ℓ−g
ℓ−1]andE[g
ℓ]vs.ℓ; (right)V[g
ℓ−g
ℓ−1]andV[g
ℓ]. QoI
g=cm(t18,x1),x1=(1.1,−0.8).
Levelℓnℓ,(
nℓ
nℓ−1
)rℓ,(
rℓ
rℓ−1
)τℓ=6016/rℓ
Computingtimes(sℓ),(
sℓ
sℓ−1
)
averagemin.max.
0 608 188 32 3 2.43.4
1 2368(3.9)376(2)16 22(7.3)15.527.8
2 9344(3.9)752(2)8 189(8.6)115237
3 37120(4)1504(2)4 1831(10)8822363
4 147968(4)3008(2)2 18580(10)786525418
Table1:#ndofsnℓ,numberoftimestepsrℓ,stepsizeintimeτℓ,average,minimal,andmaxi-
malcomp.times.
Estimatedtheweakandstrongconvergenceratesαandβ,andconstantsC1andC2fordifferentQoIs. QoIareso-
lutioninapointxiandanintegralIi(t,ω):=
R
x∈∆i
cm(t,x,ω)ρ(cm(t,x,ω))dx,∆i:=[xi−0.1,xi+0.1]×[yi−0.1,yi+0.1],
i=1,...,6.
QoI αc1β C2
cm(x1,15τ)-1.070.47-1.974·10
−3
I2(cm)-1.58.3-2.60.4
εm0m1m2m3m4
0.2280000
0.1101000
0.05574100
0.02534222410
0.01327820535610 10 20 30 40
0.05
0.1
0.15
0.2
0.25
0.3
0.35 0 10 20 30 40
0
1
2
3
4
5
10
-4 0 10 20 30 40
-8
-6
-4
-2
0
10
-3 0 10 20 30 40
0
0.2
0.4
0.6
0.8
1
10
-6
Figure5:(1)ThecoefficientofvarianceCV
ℓ:=CV(g
ℓ),(2)thevarianceV[g
ℓ],(3)themeanE[g
ℓ−g
ℓ−1],(4)the
varianceV[g
ℓ−g
ℓ−1]. The QoI isg=cm(t,x1). The small oscillations in the two lower pictures are due to the
dependenceoftherechargeˆqxonthetime,cf.(7).10
-2
10
-1
10
0
10
5
10
10
ComplexitycomparisonofMLandMLMCvs. accuracy ε(horizontalaxis)inlog-logscale.
Acknowledgements:AlexandervonHumboldtfoundationandKAUSTHPC.
1. D. Logashenko, A. Litvinenko, R. Tempone, E. Vasilyeva, G. Wittum, Uncertainty quantification in the Henry problem using the multilevel Monte Carlo method,
Journal of Computational Physics, Vol. 503, 2024, 112854,https://doi.org/10.1016/j.jcp.2024.112854
2. D. Logashenko, A. Litvinenko, R. Tempone, G. Wittum, Estimation of uncertainties in the density driven flow in fractured porous media using MLMC,
arXiv:2404.18003, 2024