Numerical Simulations of CO 2 Geo-Sequestration using PETSc Henrik - - PowerPoint PPT Presentation

SLIDE 1

Numerical Simulations of CO2 Geo-Sequestration using PETSc

Henrik B¨ using

Institute for Applied Geophysics and Geothermal Energy E.ON Energy Research Center RWTH Aachen University

June 30th, 2016

SLIDE 2

Two-phase flow Numerical method and test example Two-phase two-component flow Properties of CO2 and brine

Overview

Two-phase flow in porous media Numerical method and test example Two-phase two-component flow Properties of CO2 and brine

• H. B¨

using Numerical Simulations of CO2 Geo-Sequestration 2

SLIDE 3

Two-phase flow Numerical method and test example Two-phase two-component flow Properties of CO2 and brine

Representative elementary volume (REV)

microscale

gas phase

REV

averaging rock matrix liquid phase

Porosity: φ = Vpores

Vtotal , Saturation of phase α: Sα = Vα Vpores ,

Absolute permeability: K = kf

µ ρg .

• H. B¨

using Numerical Simulations of CO2 Geo-Sequestration 3

SLIDE 4

Two-phase flow Numerical method and test example Two-phase two-component flow Properties of CO2 and brine

Initial-Boundary-Value problem

pw-Sn-formulation ∂(φρw(1 − Sn)) ∂t + div

• ρw

krw(Sn) µw K(∇pw − ρwg)

• = ρwqw

∂(φρnSn) ∂t + div

• ρn

krn(Sn) µn K(∇pw + ∇pc(Sn) − ρng)

• = ρnqn

Initial conditions Sn(x, 0) = Sn0(x), pw(x, 0) = pw0(x) x ∈ Ω Boundary conditions pw(x, t) = gDw(x, t) on ΓDw ρwv w · n = gNw(x, t) on ΓNw Sn(x, t) = gDn(x, t) on ΓDn ρnv n · n = gNn(x, t) on ΓNn

• H. B¨

using Numerical Simulations of CO2 Geo-Sequestration 4

SLIDE 5

Two-phase flow Numerical method and test example Two-phase two-component flow Properties of CO2 and brine

Nonlinearities

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

water saturation Sw [−] relative permeability kr [−]

Brooks−Corey, k

rw

Brooks−Corey, k

rn

van Genuchten krw van Genuchten krn 0.2 0.4 0.6 0.8 1 1 2 3 4 5 6 7 8 9 10 x 10

5

water saturation Sw [−] capillary pressure pc [Pa]

Brooks−Corey van Genuchten

Brooks-Corey krw = S

2+3λ λ

e

krn = (1 − Se)2

• 1 − S

2+λ λ

e

• pc = pdS−1/λ

e

van Genuchten krw = √ Se

• 1 − (1 − S1/m

e

)m2 krn = (1 − Se)

1 3

• 1 − S

1 m

e

2m pc = 1 α(S−1/m

e

− 1)1/n

Effective saturation: Se =

Sw−Swr 1−Swr −Snr ,

0 ≤ Se ≤ 1

• H. B¨

using Numerical Simulations of CO2 Geo-Sequestration 5

SLIDE 6

Two-phase flow Numerical method and test example Two-phase two-component flow Properties of CO2 and brine

Numerical method

∂(φραSα) ∂t + div

• ρα

krα µα K(∇pα − ραg)

• = ραqα

α ∈ {w, n}

◮ First step: Semidiscretization in space with two-point flux

• approximation. Leads to a system of ordinary differential

equations.

◮ Second step: Time-Integration with implicit Euler method.

Leads to a system of nonlinear algebraic equations (remember relative permeabilities and capillary pressure). F(u) = 0 with u =  pw Sn   and F =  F1 F2   Linearize this nonlinear system of equations with Newton’s method.

• H. B¨

using Numerical Simulations of CO2 Geo-Sequestration 6

SLIDE 7

Two-phase flow Numerical method and test example Two-phase two-component flow Properties of CO2 and brine

Numerical method

• α

φ(ραSα)n+1

i

− (ραSα)n

i

∆t Vi +

• α
• j
• ρα

krα µα K n+1

ij

pw,j − pw,i di + dj − ρijgij n+1 Aij −

• α

qn+1

α,i Vi = 0

Two-point flux approximation for two neighbouring grid cells i and j with distances di and dj to the interface separating the two control volumes with area Aij.

• H. B¨

using Numerical Simulations of CO2 Geo-Sequestration 7

SLIDE 8

Two-phase flow Numerical method and test example Two-phase two-component flow Properties of CO2 and brine

Newton’s method

Transformation into linear system ∂F(u) ∂u ∆u = −F(u) Jacobian J := ∂F(u)

∂u

and ∆u := uj+1 − uj. Jacobian is of the form J =  

∂F1 ∂pw ∂F1 ∂Sn ∂F2 ∂pw ∂F2 ∂Sn

  Exact Jacobian computed by Automatic Differentiation (AD) using ADiMat, TAPENADE

• r TAF.

Every quadrant has non-zero entries due to coupling of equations.

• H. B¨

using Numerical Simulations of CO2 Geo-Sequestration 8

SLIDE 9

Two-phase flow Numerical method and test example Two-phase two-component flow Properties of CO2 and brine

Comparison of exact and approximate Jacobians

Jij = ∂Fi(u) ∂uj ≈ Fi(. . . , uj−1 + ∆uj, uj+1, . . .) − Fi(. . . , uj−1 − ∆uj, uj+1, . . .) 2∆uj with u = (pw, Sn)T = (u1, u2, . . . , uN)T and ∆uj = δ · uj.

1 10 20 30 40 50 60 70 80 90 5 10 15 20

Time step number # Newton iterations

Homogeneous case

Finite differences (FD) Automatic differentiation (AD)

Exact Jacobians save time: One vs. two evaluations. Newton iterations decrease.

• H. B¨

using Numerical Simulations of CO2 Geo-Sequestration 9

SLIDE 10

Two-phase flow Numerical method and test example Two-phase two-component flow Properties of CO2 and brine

Comparison of exact and approximate Jacobians

Jij = ∂Fi(u) ∂uj ≈ Fi(. . . , uj−1 + ∆uj, uj+1, . . .) − Fi(. . . , uj−1 − ∆uj, uj+1, . . .) 2∆uj with u = (pw, Sn)T = (u1, u2, . . . , uN)T and ∆uj = δ · uj.

1 10 20 30 40 50 60 70 80 90 5 10 15 20

Time step number # Newton iterations

Heterogeneous case

Finite differences (FD) Automatic differentiation (AD)

Exact Jacobians save time: One vs. two evaluations. Newton iterations decrease.

• H. B¨

using Numerical Simulations of CO2 Geo-Sequestration 10

SLIDE 11

Two-phase flow Numerical method and test example Two-phase two-component flow Properties of CO2 and brine

Used preconditioners and iterative solvers Balay et al. (1997)

Algebraic multigrid

◮ Hypre/BoomerAMG http://acts.nersc.gov/hypre/ ◮ Notay (2012)/AGMG

http://homepages.ulb.ac.be/~ynotay/AGMG/

◮ PETSc/GAMG http://www.mcs.anl.gov/petsc/ ◮ Trilinos/ML http://trilinos.sandia.gov/packages/ml/

Solvers

◮ MUMPS/LU

http://graal.ens-lyon.fr/MUMPS/

◮ BiCGStab ◮ GMRES ◮ FGMRES ◮ Geometric multigrid

(2 and 3 level) Preconditioners

◮ Incomplete LU ◮ Hypre/Euclid ◮ Block-Jacobi

• H. B¨

using Numerical Simulations of CO2 Geo-Sequestration 11

SLIDE 12

Two-phase flow Numerical method and test example Two-phase two-component flow Properties of CO2 and brine

Heterogeneous porosity and permeability

Gaussian distribution for Porosity field. Permeability after Pape et al. (1999). Fractal model valid for Rotliegend sandstone of NE-German basin: K = 155 φ + 37315 φ2 + 630(10 φ)10.

• H. B¨

using Numerical Simulations of CO2 Geo-Sequestration 12

SLIDE 13

Two-phase flow Numerical method and test example Two-phase two-component flow Properties of CO2 and brine

Performance of iterative solvers and preconditioners

100 110 120 130 140 150 10

−14

10

−12

10

−10

10

−8

10

−6

10

−4

Time [%] log10(Residual)

Homogeneous porous medium

MUMPS/LU Hypre/Euclid ILU0 ASM Trilinos/ML (W Cycle) Trilinos/ML (V Cycle) 100 120 140 160 180 200 220 10

−14

10

−12

10

−10

10

−8

10

−6

10

−4

Time [%] log10(Residual)

Heterogeneous porous medium

MUMPS/LU Hypre/Euclid ILU0 ASM Trilinos/ML (W Cycle) Trilinos/ML (V Cycle) 90 95 100 105 110 115 120 125 130 10

−14

10

−12

10

−10

10

−8

10

−6

10

−4

Time [%] log10(Residual)

Homogeneous porous medium

MUMPS/LU BiCGStab+Hypre/Euclid Geometric MG (2 level) Geometric MG (3 level) FGMRES+Hypre/Euclid GMRES+Hypre/Euclid 100 110 120 130 140 10

−14

10

−12

10

−10

10

−8

10

−6

10

−4

Time [%] log10(Residual)

Heterogeneous porous medium

MUMPS/LU BiCGStab+Hypre/Euclid Geometric MG (2 level) Geometric MG (3 level) FGMRES+Hypre/Euclid GMRES+Hypre/Euclid

Geometric multigrid best. Necessity for large-scale problems.

• H. B¨

using Numerical Simulations of CO2 Geo-Sequestration 13

SLIDE 14

Two-phase flow Numerical method and test example Two-phase two-component flow Properties of CO2 and brine

CO2 injection into heterogeneous porous media.

x-Extension [m] x-Extension [m]

z-Extension [m] z-Extension [m]

Day: Day:

• H. B¨

using Numerical Simulations of CO2 Geo-Sequestration 14

SLIDE 15

Two-phase flow Numerical method and test example Two-phase two-component flow Properties of CO2 and brine

Convergence study

Grid size: I0·J0·K0 = (2x·6+1)·2·(2x+1) EOCi+1 = log(2)−1|log ei ei+1

• |

x Nodes MUMPS/LU [s] ILU0 [s] GeoMG3 [s] EOC(pw) 2 250 114 75 106 1.32

• H. B¨

using Numerical Simulations of CO2 Geo-Sequestration 15

SLIDE 16

Two-phase flow Numerical method and test example Two-phase two-component flow Properties of CO2 and brine

Convergence study

Grid size: I0·J0·K0 = (2x·6+1)·2·(2x+1) EOCi+1 = log(2)−1|log ei ei+1

• |

x Nodes MUMPS/LU [s] ILU0 [s] GeoMG3 [s] EOC(pw) 2 250 114 75 106 1.32 3 882 374 403 340 0.99

• H. B¨

using Numerical Simulations of CO2 Geo-Sequestration 16

SLIDE 17

Two-phase flow Numerical method and test example Two-phase two-component flow Properties of CO2 and brine

Convergence study

Grid size: I0·J0·K0 = (2x·6+1)·2·(2x+1) EOCi+1 = log(2)−1|log ei ei+1

• |

x Nodes MUMPS/LU [s] ILU0 [s] GeoMG3 [s] EOC(pw) 2 250 114 75 106 1.32 3 882 374 403 340 0.99 4 3298 1396 1533 1262 1.00

• H. B¨

using Numerical Simulations of CO2 Geo-Sequestration 17

SLIDE 18

Two-phase flow Numerical method and test example Two-phase two-component flow Properties of CO2 and brine

Convergence study

Grid size: I0·J0·K0 = (2x·6+1)·2·(2x+1) EOCi+1 = log(2)−1|log ei ei+1

• |

x Nodes MUMPS/LU [s] ILU0 [s] GeoMG3 [s] EOC(pw) 2 250 114 75 106 1.32 3 882 374 403 340 0.99 4 3298 1396 1533 1262 1.00 5 12738 5899 7270 5339 1.00

• H. B¨

using Numerical Simulations of CO2 Geo-Sequestration 18

SLIDE 19

Two-phase flow Numerical method and test example Two-phase two-component flow Properties of CO2 and brine

Convergence study

Grid size: I0·J0·K0 = (2x·6+1)·2·(2x+1) EOCi+1 = log(2)−1|log ei ei+1

• |

x Nodes MUMPS/LU [s] ILU0 [s] GeoMG3 [s] EOC(pw) 2 250 114 75 106 1.32 3 882 374 403 340 0.99 4 3298 1396 1533 1262 1.00 5 12738 5899 7270 5339 1.00 6 50050 28350 45476 22336 1.00

• H. B¨

using Numerical Simulations of CO2 Geo-Sequestration 19

SLIDE 20

Two-phase flow Numerical method and test example Two-phase two-component flow Properties of CO2 and brine

Convergence study

Grid size: I0·J0·K0 = (2x·6+1)·2·(2x+1) EOCi+1 = log(2)−1|log ei ei+1

• |

x Nodes MUMPS/LU [s] ILU0 [s] GeoMG3 [s] EOC(pw) 2 250 114 75 106 1.32 3 882 374 403 340 0.99 4 3298 1396 1533 1262 1.00 5 12738 5899 7270 5339 1.00 6 50050 28350 45476 22336 1.00 7 198402 152108 189875 187140

• H. B¨

using Numerical Simulations of CO2 Geo-Sequestration 20

SLIDE 21

Two-phase flow Numerical method and test example Two-phase two-component flow Properties of CO2 and brine

Two-phase two-component flow

• α∈{w,n}

∂(φραxκ

αSα)

∂t −

• α∈{w,n}

div(ραλαxκ

αK(∇pα − ραg)

−

• α∈{w,n}

div(ραDκ

pm,α∇xκ α) − qκ = 0,

κ ∈ {H2O, CO2} (2p2c) Special case: Two-phase flow xCO2

n

= 1, xH2O

n

= 0 xCO2

w

= 0, xH2O

w

= 1 ∂φρwSw ∂t − div(ρwλwK(∇pw − ρwg)) = qw ∂φρnSn ∂t − div(ρnλnK(∇pn − ρng)) = qn (2p)

• H. B¨

using Numerical Simulations of CO2 Geo-Sequestration 21

SLIDE 22

Two-phase flow Numerical method and test example Two-phase two-component flow Properties of CO2 and brine

Closure relations and primary variables

Algebraic closure relations:

• α∈{w,n}

Sα = 1, pc = pn − pw

• c∈{H2O,CO2}

xc

α = 1,

α ∈ {w, n} Choose primary variables: pw, Sn. Dependent variables: xc

α = xc α(pn, T, sal), ρα = ρα(pα, T, sal, xc α), µα = µα(pα, T, sal).

• H. B¨

using Numerical Simulations of CO2 Geo-Sequestration 22

SLIDE 23

Two-phase flow Numerical method and test example Two-phase two-component flow Properties of CO2 and brine

Phase diagram

Phase diagram for two- component system

A B Sn = 0 Sn = 1 0 < Sn < 1 p Total concentration of CO2 Equilibrium single-phase liquid Equilibrium two-phase zone Equilibrium single-phase gas

xc

α gives mole of component c per total mole in phase α when the

two phases are in equilibrium. Problem: Equations only hold for two-phase regions. Not in single-phase regions. Limit of equations for Sn → 0: ∂(φρwxc

w)

∂t − div( ρw µw xc

wK(∇pw − ρwg))

− div(ρwDc

w∇xc w) + qc = 0,

c ∈ {H2O, CO2} (2c)

• H. B¨

using Numerical Simulations of CO2 Geo-Sequestration 23

SLIDE 24

Two-phase flow Numerical method and test example Two-phase two-component flow Properties of CO2 and brine

Extended Saturations

Solution:

◮ Introduce residual saturations and avoid single-phase regions

→ unrealistic.

◮ Switch primary variables, choose e.g. xCO2 w

and pw.

◮ Extend concept of saturation and use two-phase flow

equations everywhere. Method of extended saturations after Abadpour & Panfilov (2008). Idea: Introduce imaginary gas phase in zone of undersaturated liquid and imaginary liquid phase for zone of oversaturated gas. ˜ S < 0 undersaturated liquid 0 ≤ ˜ S ≤ 1 in the two-phase region ˜ S > 1

• versaturated gas

Sn =      if ˜ S < 0 ˜ S if 0 ≤ ˜ S ≤ 1 1 if ˜ S > 1.

• H. B¨

using Numerical Simulations of CO2 Geo-Sequestration 24

SLIDE 25

Two-phase flow Numerical method and test example Two-phase two-component flow Properties of CO2 and brine

Consistence conditions

Consistence conditions for imaginary gas: ˜ S < 0 (undersaturated liquid). ρn = ρw, µn = µw krw( ˜ S) = 1 − ˜ S, krn( ˜ S) = ˜ S pc( ˜ S) = 0 Dn = Dw

• 1 + xCO2

n

− xCO2

w

˜ S ∇ ˜ S ∇−1xCO2

n

• xCO2

n

= xCO2

n

(pn, T), xCO2

w

= xCO2

w

(pn, T) Plugging consistence equations into (2p2c) leads to correct single-phase equations (2c).

• H. B¨

using Numerical Simulations of CO2 Geo-Sequestration 25

SLIDE 26

Two-phase flow Numerical method and test example Two-phase two-component flow Properties of CO2 and brine

Density

Density of CO2: ρn = ρn(pn, T)

(Span & Wagner, 1996)

5 10 15 20 100 200 300 400 500 600 700 800 900 1000

Pressure [MPa] Density [kg m−3]

CO2 density

5 15 25 35 45 65 85 Temperature [°C]

Density of brine: ρw = ρw(pw, T, sal, xCO2

w

)

(Batzle & Wang, 1992; Garcia, 2001)

10 20 30 40 50 60 70 80 90 100 950 1000 1050 1100 1150 1200 1250

Temperature [°C] Density [kg m−3]

Water density

pure water brine brine with dissolved CO2

Pressure: pw = 10 MPa Salinity: sal = 0.25 mol mol−1 Dissolved CO2: xCO2

w

= 0.02 mol mol−1

• H. B¨

using Numerical Simulations of CO2 Geo-Sequestration 26

SLIDE 27

Two-phase flow Numerical method and test example Two-phase two-component flow Properties of CO2 and brine

Viscosity

Viscosity of CO2: µn = µn(pn, T)

(Fenghour et al., 1998)

5 10 15 20 20 30 40 50 60 70 80 90 100 110 120

Pressure [MPa] Viscosity [10−6 Pa s]

CO2 viscosity

5 15 25 35 45 65 85 Temperature [°C]

Viscosity of brine: µw = µw(T, sal)

(Batzle & Wang, 1992)

10 20 30 40 50 60 70 80 90 100 0.5 1 1.5 2 2.5

Temperature [°C] Viscosity [10−3 Pa s]

Brine viscosity

0.00 0.05 0.10 0.15 0.20 0.25 Salinity [mol/mol]

• H. B¨

using Numerical Simulations of CO2 Geo-Sequestration 27

SLIDE 28

Two-phase flow Numerical method and test example Two-phase two-component flow Properties of CO2 and brine

Solubility

Solubility of CO2 in brine: xCO2

w

= xCO2

w

(pn, T, sal)

(Spycher et al., 2005)

Solubility of H2O in gas: xH2O

n

= xH2O

n

(pn, T, sal)

(Spycher et al., 2005)

10 20 30 40 50 60 0.5 1 1.5 2 2.5 3 3.5

CO2 mole fraction in wetting phase

Pressure [MPa] xCO

2

w

× 100 [mol CO2 / mol H2O]

1 2 4 Salinity [mol/kg] 10 20 30 40 50 60 1 2 3 4 5 6 7 8 9 10

H2O mole fraction in non−wetting phase

Pressure [MPa] xH

2O

n

× 1000 [mol H2O / mol CO2]

1 2 4 Salinity [mol/kg]

Temperature: T = 30 ◦C, Salinity: Different molalities of NaCl.

• H. B¨

using Numerical Simulations of CO2 Geo-Sequestration 28

SLIDE 29

Two-phase flow Numerical method and test example Two-phase two-component flow Properties of CO2 and brine

Numerical simulation of CO2 injection.

• H. B¨

using Numerical Simulations of CO2 Geo-Sequestration 29

SLIDE 30

Two-phase flow Numerical method and test example Two-phase two-component flow Properties of CO2 and brine

Summary and conclusion

Summary:

◮ Test of preconditioners and iterative solvers ◮ CO2 injection into highly heterogeneous porous media ◮ Convergence study ◮ Comparison of automatic differentiation (AD) and finite

differences (FD) Conclusion:

◮ Difficulties with algebraic multigrid due to hyperbolic

character of equations

◮ Geometric multigrid performs favorable ◮ Linear increase of computation time ◮ AD outperforms FD in terms of precision and speed

Thank you for your attention!

• H. B¨

using Numerical Simulations of CO2 Geo-Sequestration 30

SLIDE 31

Two-phase flow Numerical method and test example Two-phase two-component flow Properties of CO2 and brine

Vector Form

Assuming constant density and porosity S ∂u ∂t − div(c∇u − G) = f with S =  0 −φρw φρn   , c =  ρwλw(Sn)K ρnλn(Sn)K ρnλn(Sn)K dpc(Sn)

dSn

  f = ρwqw ρnqn

• ,

G = ρwλw(Sn)Kρwg ρnλn(Sn)Kρng

• and u =

pw Sn

• .
• H. B¨

using Numerical Simulations of CO2 Geo-Sequestration 31

SLIDE 32

Two-phase flow Numerical method and test example Two-phase two-component flow Properties of CO2 and brine

References I

Balay, S., Gropp, W. D., McInnes, L. C., & Smith, B. F., 1997. Efficient management of parallelism in object oriented numerical software libraries, in E. Arge, A. M. Bruaset, & H. P. Langtangen (eds.), Modern Software Tools in Scientific Computing, pp. 163–202, Birkh¨ auser Press. Batzle, M. & Wang, Z., 1992. Seismic properties of pore fluids, Geophysics, 57(11), 1396–1408. Brooks, R. J. & Corey, A. T., 1964. Hydraulic properties of porous media, vol. 3, Colorado State University Hydrology Paper, Fort Collins. Fenghour, A., Wakeham, W. A., & Vesovic, V., 1998. The viscosity of carbon dioxide, Journal of Physical and Chemical Reference Data, 27(1), 31–44.

• H. B¨

using Numerical Simulations of CO2 Geo-Sequestration 32

SLIDE 33

Two-phase flow Numerical method and test example Two-phase two-component flow Properties of CO2 and brine

References II

Garcia, J. E., 2001. Density of aqueous solutions of CO2, Tech. rep., Earth Sciences Division, Lawrence Berkeley National Laboratory. Griewank, A., 2000. Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation, Society of Industrial and Applied Mathematics (SIAM), Philadelphia, PA. Notay, Y., 2012. Aggregation-based algebraic multigrid for convection-diffusion equations, SIAM Journal on Scientific Computing, 34, A2288–A2316. Pape, H., Clauser, C., & Iffland, J., 1999. Permeability prediction based on fractal pore-space geometry, Geophysics, 64(5), 1447–1460.

• H. B¨

using Numerical Simulations of CO2 Geo-Sequestration 33

SLIDE 34

Two-phase flow Numerical method and test example Two-phase two-component flow Properties of CO2 and brine

References III

Span, R. & Wagner, W., 1996. A new equation of state for carbon dioxide covering the fluid region from the triple-point temperature to 1100 K at pressures up to 800 MPa, Journal of Physical and Chemical Reference Data, 25(6), 1509–1596. Spycher, N., Pruess, K., & Ennis-King, J., 2005. CO2-H2O mixtures in the geological sequestration of CO2. II. partitioning in chloride brines at 12–100 ◦C and up to 600 bar, Geochimica et Cosmochimica Acta, 69(13), 3309–3320. van Genuchten, M. T., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils, Soil Science Society of America, 44, 892–898.

• H. B¨

using Numerical Simulations of CO2 Geo-Sequestration 34