Heft 208 Nguyen Nghia Hung
Sediment dynamics in the floodplain
of the Mekong Delta, Vietnam
Institut für Wasser- und Umweltsystemmodellierung
Heft 212 Alexandru-Bogdan Tatomir
From Discrete to Continuum Concepts
in Fractured Porous Media
Von der Fakultät Bau- und Umweltingenieurwissenschaften der
Universität Stuttgart zur Erlangung der Würde eines
Doktor-Ingenieurs (Dr.-Ing.) genehmigte Abhandlung
Vorgelegt von
Alexandru-Bogdan Tatomir
aus Bukarest, Rumänien
Hauptberichter: Prof. Dr.-Ing. Rainer Helmig Mitberichter: Prof. Dr. Sorin Pop
Habil. Prof. Dr.-Ing. Holger Class Asst. Prof. Adam Szymkiewicz Tag der mündlichen Prüfung: 15. Februar 2012
Institut für Wasser- und Umweltsystemmodellierung
der Universität Stuttgart
Heft 212
From Discrete to Continuum
Concepts of Flow in Fractured
Porous Media
von
Dr.-Ing.
Alexandru-Bogdan Tatomir
Eigenverlag des Instituts für Wasser- und Umweltsystemmodellierung
der Universität Stuttgart
Bibliografische Information der Deutschen Nationalbibliothek
Die Deutsche Nationalbibliothek verzeichnet diese Publikation in der Deutschen Nationalbibliografie; detaillierte bibliografische Daten sind im Internet über http://www.d-nb.de abrufbar
Tatomir, Alexandru-Bogdan:
From Discrete to Continuum Concepts of Flow in Fractured Porous Media / von Alexandru-Bogdan Tatomir. Institut für Wasser- und
Umweltsystemmodellierung, Universität Stuttgart. - Stuttgart: Institut für Wasser- und Umweltsystemmodellierung, 2012
(Mitteilungen Institut für Wasser- und Umweltsystemmodellierung, Universität Stuttgart: H. 212)
Zugl.: Stuttgart, Univ., Diss., 2012 ISBN 978-3-942036-16-0
NE: Institut für Wasser- und Umweltsystemmodellierung <Stuttgart>: Mitteilungen Gegen Vervielfältigung und Übersetzung bestehen keine Einwände, es wird lediglich um Quellenangabe gebeten.
Herausgegeben 2012 vom Eigenverlag des Instituts für Wasser- und Umwelt-systemmodellierung
This thesis is dedicated to:
My beloved mother who as a physicist inspired me with her love for science. She will always remain in my heart.
My beloved father who always supported and encouraged me on this way. My two little daughters for their unending resource of joy.
My adorable wife for her constant support, caring and patience.
My family for their emotional support and for always being besides me when I most need them.
I would like to thank all those who encouraged and helped me to pursue a PhD. First of all, I want to thank my main supervisor Prof. Rainer Helmig for the opportunity to become a member of his esteemed group, for his creative ideas and interesting discussions, for his engagement, positive attitude and fairness. I would also like to express my gratitude to my co-supervisor Prof. Sorin Pop for his constant support and valuable advices. I want to thank Dr. Holger Class for carefully reviewing my dissertation and for the advices. I want to express my deepest appreciation to Dr. Adam Szymkiewicz for his help, guidance and friendship.
I am also very grateful to the doctoral programme “Environment and Water” (ENWAT) and the scholarship program IPSWaT.
I want to thank my colleagues at the IWS, particularly LH2, for the scientific and non-scientific help and the pleasant working environment.
I convey my appreciation to Bernd Flemisch and Andreas Lauser for their support to establish the models.
The period in the LH2 group was very enjoyable also due to the football groups of IWS, IAHR and ENWAT.
I would like to express my deepest appreciation to Dieter and Gina for their tremendous help, my friends George, Veronica and Mihaela with whom we spent enjoyable weekends together.
I thank the ”Friday lunch group at the Chinese”, Onur, Andreas, Jeff and Hab-tamu, which brought spice to our PhD time.
Moreover, I would like to thank Prof. Mirel and my former colleagues from the Technical University of Civil Engineering Bucharest, especially Prof. Sandu Lucian and Prof. Ioan Bica for their support.
Nomenclature XII
Abstract XIV
Zusammenfassung XVI
1 Introduction 1
1.1 Motivation . . . 1
1.2 Objectives of the Research . . . 4
1.3 Thesis Outline . . . 5
1.4 Structure . . . 6
2 Models of Fluid Flow in Fractured Media 8 2.1 Scales of Consideration . . . 8
2.2 Properties of Porous Media . . . 9
2.2.1 Fundamental Definitions . . . 9 2.2.2 Fluid Properties . . . 11 2.2.2.1 Density . . . 11 2.2.2.2 Viscosity . . . 11 2.2.3 Matrix Properties . . . 12 2.2.3.1 Porosity . . . 12 2.2.3.2 Intrinsic Permeability . . . 12
2.2.4 Fluid-Matrix Interaction Properties . . . 13
2.2.4.1 Saturation and Residual Saturation . . . 13
2.2.4.2 Interfacial Tension and Wettability . . . 14
2.2.4.3 Capillary Pressure . . . 14
2.2.4.4 Relative Permeability . . . 16
2.3 General Equations of Two-Phase Fluid Flow in Porous Media . . . . 18
2.3.1 The Darcy Equation . . . 19
2.3.2 Mass Conservation Equation . . . 19
2.3.3 Phase Pressure-Saturation Formulation . . . 21
2.4 Overview of Fractured Systems . . . 22
2.4.1 Scales of Consideration for Fractured Media . . . 22
2.4.3 Fracture Geometry Analysis . . . 25
2.4.4 Geometrical Models for Single Open Fractures . . . 28
2.5 Discrete Fracture Models (DFM) . . . 30
2.5.1 Discrete Model Approaches . . . 30
2.5.2 Lower and Equi-Dimensional Approach . . . 31
2.5.3 Relative Permeability - Saturation Relation in Fractures . . . . 32
2.5.4 Capillary Pressure in Fractures . . . 34
2.5.5 Interface Conditions at Media Discontinuities . . . 35
2.5.6 Mathematical Model of the DFM . . . 37
2.5.6.1 Lower-dimensional approach . . . 39
2.5.6.2 Equi-dimensional approach . . . 41
2.6 Continuum Fracture Models (CFM) . . . 43
2.6.1 Overview . . . 43
2.6.2 Matrix-Fracture Exchange . . . 45
2.6.3 Limitations of the Standard Double-Porosity Models . . . 46
2.7 Multiple Interacting Continua (MINC) Method . . . 47
2.7.1 General Overview of the Method . . . 47
2.7.2 The Geometric Description of the MINC Block . . . 49
2.7.3 Mathematical Model for the MINC . . . 50
2.7.4 Further Developments of MINC Method and Additional Con-siderations . . . 50
2.8 Extended MINC Method . . . 52
2.8.1 Matrix-block Subgridding - Intra-coarse-block Transmissibilities 52 2.8.2 Coarse Blocks Inter-transmissibilities . . . 55
3 Numerical Implementation 57 3.1 Numerical Modeling . . . 57
3.1.1 Spatial and Temporal Discretization . . . 57
3.1.1.1 Dual Meshing . . . 57
3.1.1.2 Weak Formulation . . . 58
3.1.1.3 Two-phase Flow Equations Discretization . . . 62
3.1.1.4 Two-phase Flow Equations Discretization for DFM . 63 3.1.1.5 Two-phase Flow Equations Discretization for Ex-tended MINC Method . . . 64
3.1.2 Linearization of the Partial Differential Equation System . . . 65
3.2 DuMuxModeling System . . . . 66
3.3 2pDFM simulator . . . 67
3.3.1 Geometry Modeling . . . 67
3.3.2 DFM Flow Model . . . 71
3.4 2pMINC simulator . . . 71
3.5 Workflow . . . 74
3.5.1 Preprocessing . . . 74
3.5.2 Processing . . . 74
3.5.3 Postprocessing . . . 76
4 Attempts at Model Verification and Validation 77 4.1 2pDFM Simulator Verification . . . 77
4.1.1 Example DFM 1: Kazemi laboratory experiment . . . 77
4.1.2 Example DFM 2: Flow parallel to a fracture: Equi- vs lower-dimensional approach (2pDFM-L vs. 2pDFM-E) . . . 82
4.1.2.1 Example DFM 2a: Vertical gas infiltration experiment into a single fully saturated fracture . . . 82
4.1.2.2 Example DFM 2b: Vertical gas infiltration experiment into a single partially saturated fracture . . . 85
4.1.2.3 Example DFM 2c: Vertical gas infiltration experiment into a single partially saturated fracture . . . 90
4.1.2.4 Example DFM 2d: Vertical water imbibition into par-tially air saturated medium . . . 91
4.1.3 Example DFM 3: Flow perpendicular to a fracture: Equi- vs lower- dimensional approach (2pDFM-L vs. 2pDFM-E) . . . . 98
4.1.4 Example DFM 4: Steady-state in a rock mass intersected by fractured zones: HYDROCOIN Case 2 Level 1 . . . 98
4.2 2pMINC Model Verification . . . 109
4.2.1 Example MINC 1: Single MINC block vs. quasi-1D domain problem . . . 109
4.2.2 Example MINC 2: Single homogeneous MINC block . . . 110
4.2.3 Example MINC 3: Single MINC block with interior heteroge-neous lens . . . 117
4.2.4 Example MINC 4: Two phase MINC model verification with an analogy to a 2D heterogeneous domain . . . 120
4.2.5 Example MINC 5: Series of square-shaped nested blocks . . . 125
4.2.6 Example MINC 6: Crossing fractures . . . 130
5 Model Application to Field Scale Problems 135 5.1 Two phase flow in an idealized periodic fracture system . . . 135
5.1.1 “Permeameter type” experiment . . . 135
5.1.2 Quarter five spot problem . . . 141
5.2 Bristol Fractured Aquifer . . . 142
5.2.1 Air-Water System . . . 144
5.2.2 CO2- Water System . . . 146 6 Summary and Conclusions 153
1 1) Gekl ¨uftetes Reservoir 2) schematische Repr¨asentation des MINC Modells (3D) 3) 2D Querschnitt durch das ”W ¨urfelzuckermodell”4) Konnektivit¨aten: a) Konnektivit¨at der grobmaßst¨ablichen Bl¨ocke oder Kluft-Kluft; b) Kluft-Matrix c) Matrix-Matrix. . . XXI 2 Darstellung des Druckprofils und der Isodrucklinien in
quasi-station¨are Zustand f ¨ur den lokalen Block Ωiauf der Feinskala und
Konstruktion der idealisierten verschachtelten Volumenelemente Ωi
1,Ωi2,...,Ωi6der W ¨urfelzuckerbl¨ocke . . . XXIV
3 Bestimmung der grobmaßst¨ablichen Zwischenblocktransmissibilit¨at der grobmaßst¨ablicher Bl¨ocke Ωiund Ωjmit mittlerem Blockdruck pi
und pj . . . XXV
4 S¨attigungsprofile der nicht-benetzenden Phase auf der Diagonalen f ¨ur das permeameter¨ahnliche Experiment. Obere Werte beziehen sich auf DFM L¨osungen, untere Werte wurden mittels des MINC Modells berechnet. . . XXVII 1.1 Examples of fracture networks in various materials . . . 2 1.2 Examples of fracture networks in various materials . . . 3 2.1 Definition of a representative elementary volume (REV) for porous
medium and a fractured system (modified after Bear et al. [1993]) . . 10 2.2 Capillary tube connecting two immiscible fluids with interfacial
ten-sions between solid and non-wetting fluid σsn, between solid and
wetting fluid σsw, between wetting and non-wetting fluid σwn. At
equilibrium the capillary forces Fcap are equal to the gravitational
forces Fg, the fluid height is h and the contact angle θ. . . 15
2.3 Capillary pressure-saturation curves fitted to the experimental mea-surements of Kazemi [1976] with Brooks-Corey and Van Genuchten models 17 2.4 Relative permeability - saturation relationships expressed with
Brooks-Corey and Van Genuchten models . . . 18 2.5 Map of fractures in Carboniferous Limestone at Sheshymore, The
Burren, Co.Claire, Ireland. Coordinates in meters (from Odling [1999]). 23 2.6 Sketch of the relation between model concepts and scales of the
2.7 a) Definition of a fracture aperture and b) idealization of the fracture aperture (from Reichenberger et al. [2004]) . . . 27 2.8 Fracture connectivity (after Singhal and Gupta [1999]). Different types
of terminations: B - blind; C - crossing ; D - diffusely connected. . . . 27 2.9 Description of the fracture modes. . . 27 2.10 Example of single fracture models (from Reichenberger et al. [2004]) 29 2.11 Comparison of relative permeability data for rough-walled open
frac-tures (from Diomampo [2001]) . . . 33 2.12 Velocity distribution of a laminar flow between two parallel plates . 34 2.13 Conceptual model of partially saturated fractured porous medium
where fractures act as barrier for water flow (Wang and Narasimhan [1985]) . . . 36 2.14 Representation of the interfaces in the schematic column experiment
(Reichenberger et al. [2003]) . . . 36 2.15 Extended capillary pressure - saturation interface condition between
domain G1 (coarse) and G2 (fine) . . . 38 2.16 Continuity of the capillary pressure leads to a discontinuity of
satura-tion across the interface . . . 38 2.17 a) Equi- and b) lower- dimensional representation of a fracture in the
DFM spatially discretized with the box method. For the DFM-L the fracture node comprises three distinct values: saturation in the upper matrix block SMup
n , saturation in the fracture SFnand saturation in the
lower matrix block SMdown
n . . . 42
2.18 Schematic diagrams of connectivity for a) porosity single-permeability model; b) the double-porosity single-single-permeability model; c) the double-porosity, double-permeability model; d) and the MINC model . . . 44 2.19 a) Naturally fractured reservoir with b) the schematic representation
of the MINC model (3D) as an idealized “sugar cube” and d) a 2D crossection showing the nested volume elements and c) the connectiv-ities between each of the continua: 1)coarse-block inter-connectivity or fracture-fracture; 2)fracture-matrix 3) matrix-matrix connections. . 48 2.20 Schematic of connectivity: a) coarse-block inter-connectivity; b)
fracture-matrix and c) matrix-matrix connections . . . 53 2.21 Pressure profile and iso-pressure lines at quasi-steady state for local
block Ωirepresented at fine-scale and the construction of the idealized
sugar-cube block’s nested volume elements Ωi
1,Ωi2,...,Ωi6 . . . 54
2.22 Inter-coarse-block transmissibility Ti jdetermination between coarse
3.1 Construction of a control volume with the BOX-method: a) Control volume Biaround node iis surrounded by set Ni= { j1, j2,.., j5} ; b)
Sub-control volume (SCV) bE1
i in element E1has barycenter G1and
the midpoints of the edge i j1is Mi j1and of the edge i j2is Mi j2. The
SCV face fE1
i j1is the segment G1Mi j1and contains the integration point
xE1
i j1where the normal nE1i j1is applied. . . 59
3.2 Modular design of DuMux(Flemisch et al. [2011b]) . . . . 66
3.3 Example of ART format and mesh generation . . . 69 3.4 Exemplification for fine and coarse scale grid construction for the
whole fractured reservoir. Fracture pattern obtained with the geo-statistical fracture generator FRAC3D . . . 70 3.5 Different approaches for determination of the nested volume elements. 72 3.6 The simulation environment . . . 74 3.7 Schematic description of the modeling system (Tatomir et al. [2011]) 75 4.1 Kazemi experiment: Experiment setup and boundary conditions . . 79 4.2 Kazemi experiment: Capillary pressure curves for fracture and matrix
for matching simulation and laboratory results (from Wu et al. [2004b]) 79 4.3 Kazemi experiment: Water imbibition into oil saturated matrix rock
simulated with 2pDFM . . . 81 4.4 Comparison of the 2pDFM model (left) with TOUGH2 simulator and
experimental data of Kazemi [1976] (right) . . . 81 4.5 Example DFM 2a: Domain geometry, boundary and initial conditions. 83 4.6 Example DFM 2a: Comparisson between the 1D fracture and 2D
fracture approaches for an example of gas infiltration in a fully water saturated media simulated with 2p and 2pDFM. . . 84 4.7 Example DFM 2a: Non-wetting flux over the cross-sections y = 0.4m
and y = 0.8m. . . 86 4.8 Example DFM 2a: Spatial non-wetting phase saturation distribution
3400 seconds for a grid refinement of 0.05m and b) a local grid refine-ment along the vertical fracture. . . 86 4.9 Domain geometry, boundary and initial conditions for Example a)
DFM 2b b) b) DFM 2c c) DFM 2d . . . 87 4.10 Example DFM 2b: Domain description for water-gas flow example in a
1 × 7 m domain. a) Adaptive gridding for the DFM-L model obtained with ART3D; b) Zoom in near the vertical fracture; c) Zoomed cross section with the fracture discretization in the DFM-E model. . . 89 4.11 Example DFM 2b (“Vertical gas injection in 90% saturated domain
with a fracture width 0.01 m”) : Comparison of saturation distribution along the vertical line x=0.5 m (through the center of the fracture) at different times obtained with lower- and equi-dimensional fracture approaches. . . 90
4.12 Example DFM 2b (“Vertical gas injection in 90% saturated domain with a fracture width 0.05 m”) : Comparison of saturation distribution along the vertical line x=0.5 m at different times obtained with lower-and equi-dimensional fracture approaches. . . 92 4.13 Example DFM 2b: Sensitivity analysis on the grid refinement.
Non-wetting phase saturation profiles along the line x = 0.5. . . 93 4.14 Example DFM 2b: Sensitivity analysis on the grid refinement. . . 94 4.15 Example DFM 2c (“Vertical gas injection in 80% saturated domain
with a fracture width 0.01 m”): Comparison of saturation distribution along the vertical line x=0.5 m at different times obtained with lower-and equi-dimensional fracture approaches. . . 95 4.16 Example DFM 2d (“Vertical water imbibition into 90% air saturated
domain with a fracture width 0.01 m”): Comparison of saturation distribution along the vertical line x=0.5 m at different times obtained with lower- and equi-dimensional fracture approaches. . . 96 4.17 Example DFM 2d: Non-wetting flux across section y = 2.8 [m] with
lower- and equi- dimensional fracture representation . . . 97 4.18 Example DFM 2d: Non-wetting phase saturation spatial distributions
after 4.0e+04 seconds with 1D (left) and 2D (right) fracture approaches 97 4.19 Example DFM 3 (Vertical gas flow perpendicular to the fracture):
spa-tial distribution of the non-wetting phase saturation after 260 seconds. . . . 99 4.20 Example DFM 3: Vertical gas flow perpendicular to the fracture:
non-wetting flux across line y = 0.6 m. . . 99 4.21 Example DFM 3: Comparison of saturation distribution along the
vertical line x = 0.5 at different times obtained with lower- and equi-dimensional fracture approaches. . . 103 4.22 Schematic description of the domain HYDROCOIN Case 2 Level 1
(from Loefman et al. [2007]) with the vertex indexes given in Table 4.6 104 4.23 Finite element grids used for the problem Hydrocoin Level 1 Case 2. 104 4.24 Spatial distribution of the hydraulic head calculated with lower- and
equi-dimensional approaches. . . 105 4.25 Example DFM 4: Comparison of the computed hydraulic heads along
horizontal lines. The results computed by DuMux and MUFTE-UG
Tatomir [2007] are presented in subfigures (a), (c), (e) while the ones computed by FEFTRA and Grundfelt [1984] (b), (d), (f). The results of the HYDROCOIN groups HYDROCOIN [1988] are presented in Figure 4.26 (c) and (d) (continuation) . . . 106 4.26 (continuation from Figure 4.25) . . . 107 4.27 Example DFM 4: Modeled path lines (stream lines) calculated with (a)
4.28 Example DFM 4: Comparison between the modeled stream lines passing through point (x=100 m, y = 800 m) with (a) 2pDFM-DuMux,
MUFTE-UG, FEFTRA simulators and (b) the HYDROCOIN groups. . . . 108 4.29 Example MINC 1: Problem setup for the verification of one MINC
block with a quasi-1D domain. T#,#+1transmissibility between
con-tinua. . . 109 4.30 Example MINC 1 (Single homogeneous MINC block vs. quasi 1D
domain): Non-wetting phase saturation Snfronts at different times
obtained with 2p and 2pMINC and the plot of Snover the line (0.0,
0.5)(10.0, 0.5) . . . 111 4.31 Example MINC 2 (Single homogeneous MINC block): Schematic
rep-resentation of the domain and the boundary conditions. . . 112 4.32 Example MINC 2: Non-wetting phase saturation profile for MINC
model discretized into 10 interacting continua (horizontal lines), the discrete 2p model (curved line) and the resulting MINC saturation (thick red line) according to the dimensions of the nested volume elements (time = 375 sec). . . 114 4.33 Example MINC 2 (Single homogeneous MINC block): Wetting
pres-sure pwprofiles at different times. . . 115
4.34 Example MINC 2 (Single homogeneous MINC block): Non-wetting phase saturation Snfronts at different times obtained with 2p (upper
left figures) and 2pMINC (lower left figures) and the plot of Snover the
diagonal line. . . 116 4.35 Example MINC 3: Schematic representation of a fine lens surrounded
by a coarse material and the equivalent representation with a MINC model . . . 117 4.36 Example MINC 3 (Single heterogeneous MINC block): Non-wetting
phase saturation Snfronts at different times obtained with 2p and
2pMINC and the plot of Snover the line (0.0, 0.5)(1.0, 0.5) . . . 119
4.37 Example MINC 4: Analogy to a 2D heterogeneous domain experiment: Schematic representation of the domain and boundary conditions . . 120 4.38 Example MINC 4: Saturation profiles for 2pMINC (upper subfigure)
and 2p (lower subfigure) simulators at different times and their plot over lines y = 0, y = 1, y = 2, y = 3 and the corresponding MINC continua123 4.39 Example MINC 4: Saturation profiles for 2pMINC (upper subfigure)
and 2p(lower subfigure) simulators at different times and their plot over lines y = 0,y = 1,y = 2, y = 3 and the corresponding MINC continua124 4.40 Example MINC 5: Problem setup for the verification MINC method
with a set of 20 serial nested blocks; structure of the porous medium for the numerical example . . . 126
4.41 Example MINC 5 (Series of square-shaped nested blocks): Non-wetting phase saturation Snfronts at different times obtained with
the explicit 2p model and the approximate 2pMINC model and the plot of Snover the diagonal line. . . 129
4.42 Example MINC 6 (Crossing fractures): problem setup . . . 130 4.43 Example MINC 6: Spatial distribution of the non-wetting phase
satura-tion Snat different times obtained with 2pDFM-L (lower-dimensional),
2pDFM-E (equi-dimensional) model approaches and 2pMINC. Snplot
over the diagonal line. . . 132 4.43 Example MINC 6: Spatial distribution of the non-wetting phase
satura-tion Snat different times obtained with 2pDFM-L (lower-dimensional),
2pDFM-E (equi-dimensional) model approaches and 2pMINC. Snplot
over the diagonal line. . . 133 4.43 Example MINC 6: Spatial distribution of the non-wetting phase
satu-ration Snat different times obtained with 2pDFM-L (lower dimensional),
2pDFM-E (equidimensional) model approaches and 2pMINC. Snplot over
the diagonal line. . . 134 5.1 Domain description and sketch of the designed experiments for
ideal-ized fracture pattern model . . . 136 5.2 Fine-scale model discretization . . . 138 5.3 a)Pressure profile and iso-pressure lines at quasi-steady state for local
block Vi; b) Vertex centered finite volume discretization of the coarse
scale model and c) The coarse-block transmissibility determination between coarse block Viand Vj . . . 139
5.4 Non-wetting phase saturation profiles plotted over the diagonal for the permeameter type experiment. Upper values the reference DFM solution, lower values obtained with MINC model . . . 140 5.5 Quarter five spot example: non-wetting phase saturation profiles
plotted over the diagonal line (white color). DFM is plotted above the MINC. . . 143 5.6 Bristol problem: Domain description and an example for the
subdivi-sion of the Bristol fracture system in coarse blocks (two cells) and their unstructured gridding. . . 146 5.7 Bristol problem: Statistical distribution of the matrix block following a
log-normal distribution (from Geiger et al. [2011]). . . 147 5.8 Bristol problem: Pressure distribution after 1e5 seconds obtained from
the intra-coarse block transmissibility determination. . . 147 5.9 Bristol air-water problem: Non-wetting phase saturation profiles
plot-ted over the diagonal for Bristol Channel fractured reservoir analogue 148 5.10 Bristol air-water problem: Non-wetting phase saturation spatial
5.11 Bristol air-water problem: Mass fluxes of wetting and non-wetting phase plotted over time at line (x = 1.8m) obtained with 2pDFM and 2pMINCsimulators . . . 150 5.12 Bristol CO2- Water system problem: Non-wetting phase saturation
profiles plotted over the diagonal for Bristol Channel fractured reser-voir analogue . . . 152
2.1 Geometrical characterization of fractures (after Reichenberger et al. [2004] according to Singhal and Gupta [1999]) . . . 26 4.1 Domain and fluid properties for Example DFM 1 (after Wu et al. [2004b]). 80 4.2 Domain and fluid properties for Example DFM 2a . . . 83 4.3 Example DFM 2a: Simulation times . . . 85 4.4 Example DFM 2b: Domain and fluid properties for “vertical nitrogen
infiltration in partially saturated 1 × 7 meter domain” . . . 88 4.5 Example DFM 2b: Grid specifications and simulation times . . . 89 4.6 Hydrocoin Case 2 Level 1: Coordinates for 1D fracture model . . . . 101 4.7 Example DFM 4 (“HYDROCOIN” Case 2 of Level 1): Domain and
fluid properties. . . 102 4.8 Domain and fluid properties for “single MINC homogeneous block”
example. . . 113 4.9 Example MINC 2 (Single MINC homogeneous block): Geometric
characteristics for the MINC model. NC - number of continuum, A#
interface area of continuum # and d#,#+1- the virtual distance from the
nearest boundary to the nested continuum # . . . 114 4.10 Example MINC 3 (Single MINC heterogeneous block): Domain and
fluid properties for. . . 118 4.11 Example MINC 4 (Analogy to a 2D heterogeneous domain): Domain
and fluid properties. . . 122 4.12 Example MINC 4: MINC parameters . . . 122 4.13 Example MINC 5 (Series of square-shaped nested blocks): Domain
and fluid properties. . . 127 4.14 Example MINC 6 (Crossing fractures): Domain and fluid properties. 131 5.1 Domain and fluid properties for “Permeameter type” problem. . . . 137 5.2 Domain and fluid properties for Quarter five-spot problem . . . 142 5.3 Domain and fluid properties for Air-Water system problem . . . 145 5.4 Simulation times and grid specifications . . . 148 5.5 Domain and fluid properties for CO2-Water system problem . . . 151
The following table shows the significant symbols used in this work. Local notations are explained in the text.
Symbol
Definition
Unit
Latin Letters:
g gravitational acceleration [m2/s]
krα relative permeability of phase α [-]
K intrinsic permeability [m2]
p pressure [Pa]
pc capillary pressure [Pa]
pd entry pressure [Pa]
qα source or sink term of the phase α [m3/(m3s) ]
Sα saturation of phase α [-]
Sαr residual saturation of phase α [-]
Se effective saturation [-]
t time [s]
Ti j coarse scale transmissibility between node i
and j [m
3]
Ti
k,k+1 inner block i transmissibility between
contin-uum k and k + 1
Greek Letters:
α phase n or w
λ empirical constant from the Brooks-Corey pc
relationship related to the pore size distribu-tion
λα phase α mobility
µ dynamic viscosity [Pa s]
∇FF coarse-scale operator (interactions between the
fracture blocks solved at global scale) ∇MM fine-scale operator (interactions inside the
coarse block between the matrix continua solved at the local scale)
φ total porosity [-]
ρ fluid density [kg/m3]
Subscripts / Superscripts:
F fracture continuum
Fi fracture continuum at the block i
i index
j nodal index
max maximum
Mk matrix continuum k
Mi
k matrix continuum k at the block i
n non-wetting phase related quantity up upwind related quantity
w wetting phase related quantity x Cartesian coordinate in x-axis y Cartesian coordinate in y-axis
As more and more engineering applications require the correct simulation of flow and transport processes in porous media, and while many of these media present a certain degree of fracturing, this work deals with the development of numerical models that can simulate two-phase flow in large-scale fractured reservoirs. Among the applications which these models are addressing to, there are the estimation of contaminant spreading and removal, the reservoir exploitation, or more recently the CO2sequestration, the geothermal reservoir exploitation, and the nuclear waste
repositories. Fractured systems are ubiquitous around the world and can occur on a variety of lengths and scales which makes difficult the development of a general model that can handle easily all of them.
The overall purpose of this work is to improve the understanding over the concepts of the multiphase flow and processes in the fractured porous media and develop a conceptual model that allows the study of two-phase flow in fractures of arbitrary size, orientation and shape.
The flow models have been roughly classified in discrete fracture models (DFM), continuum fracture models (CFM) and hybrid. The DFMs consider the fractures explicitly and therefore require huge computation power, whereas CFMs require the determination of a representative elementary volume (REV), the appropriate effective parameters and transfer functions between continua. For a large scale problem, like a CO2storage reservoir, there can be millions of fractures which might
have to be considered and could be a formidable task for a DFM simulator. In this sense a continuum model, which is in this case in the form of a generalized dual-porosity representation, does not require the fine discretization of a DFM and the detailed fracture characterization during simulations. Thus, another goal of this work is to build a reliable large scale multiphase flow simulator based on the continuum approach. Two flow simulators, 2pDFM and 2pMINC, have been developed and tested based on the two different fracture model concepts. Both simulators are implemented in the numerical toolbox DuMux. The 2pDFM model simulates the
two-phase flow in fractured porous media using a DFM approach, with a lower dimensional representation for the fracture network (DFM-L). The model is capable to account for storage in the fractures. The 2pMINC model simulates two-phase flow in fractured porous media using the multiple interacting continua (MINC) method with an improved upscaling technique. The complex transient behavior of the flow
processes in fractured porous media is captured by subgridding the coarse blocks in nested volume elements which have effective properties calculated from the detailed representation of the fracture system. In this way, the physically based approach is kept, the accuracy of the model is preserved , the common use of empirically derived transfer functions is avoided and the complexity of the problem is considerably reduced which is reflected in the speedup factors up to 1000. This research extends the applicability range of the upscaling procedure to include entry pressure effects. Moreover, a general workflow has been developed for the numerical simulation of the two-phase flow in large-scale fractured porous media.
The results are verified, validated and tested in a fully comprehensive way for both models. To test the behaviour of the simulator for field scale problems they are applied to an idealized medium with periodic fracture pattern and to a real, naturally fractured reservoir from the Bristol Channel. The evaluation shows that the extended MINC model is able to reproduce both, the large-scale permeability and the dynamics of the fracture-matrix mass transfer, correctly.
The extended MINC method and the simulation procedure is flexible as it allows choosing the accuracy of the solution, the computation speed, and allows working with spatial information about the fracture system of various complexity and detail.
Einleitung
Kluftsysteme sind auf dieser Welt allgegenw¨artig und werden seit mehr als 60 Jahren in unterschiedlichen Bereichen, wie der Hydrologie, der Erd¨olf¨orderung oder der geothermalen Energie untersucht.
F ¨ur eine Vielzahl von Anwendungsbereichen sind pr¨adiktive Simulationen der Mehr-phasenstr¨omung und des Transports in gekl ¨ufteten por¨osen Medien von großem Interesse, z.B. bei der Erschließung von Wasserspeichern zur Wasserversorgung, bei der Sanierung wasserf ¨uhrender Schichten (Aquifere) (Berkowitz [2002],Niessner et al. [2005]), bei der Erschließung von Erd¨olspeichern (Lemonnier and Bourbiaux [2010], Kazemi [1976]), bei der Erschließung geothermaler Speicher, der W¨armespeicherung, bei Prozessen in Bergbau und Mineralisierung (lokale Extraktion und Bestimmung von Erzvorkommen), bei der CO2Speicherung in geologischen Schichten (Carneiro
[2009], Kopp [2009]), bei der Atomm ¨ulllagerung (Bodvarsson et al. [1999], Reichen-berger et al. [2003]) usw.. All diese Anwendungen k¨onnen mittels Gestaltung und
¨Uberwachung durch Computermodelle besser verstanden werden.
Rechnergest ¨utzte Modellierung ist der Prozess der Entwicklung mathematischer Repr¨asentationen einer echten oder hypothetischen Situation, um diese besser zu untersuchen, und um einen Einblick in die Funktionsweise und Handhabung dieser Systeme zu erhalten. Im Falle der Endlagerung radioaktiven M ¨ulls m ¨ussen die unterirdischen Anlagen eine Sicherheitsgarantie von hunderten oder gar tausenden von Jahren bieten. Daher liegt der Schl ¨ussel zum Verst¨andnis mehrphasiger Prozesse und zur Vorhersage ¨uber das Verhalten des Reservoirs in der Erstellung eines genauen Modells. Einige der gr¨oßten Herausforderungen bei der Erstellung von Mehrphasen-Str¨omungsmodellen in gekl ¨ufteten por¨osen Medien stellen die großen Unterschiede der Klufteigenschaften und der Eigenschaften des umgebenden Gesteins, sowie die unterschiedlichen Zeitmaßst¨abe der Str¨omungsprozesse dar.
In dieser Arbeit werden mathematische Modelle und numerische Simulationswerk-zeuge f ¨ur gekl ¨uftete por¨ose Medien entwickelt, die eine wichtige Rolle bei der Un-tersuchung unterirdischer Reservoirs, bei der Gestaltung und UnUn-tersuchung von Feldversuchen und bei der Verbesserung und Optimierung der Handhabung geolo-gischer Reservoire spielen.
Zielstellung der Arbeit
Die prim¨aren Ziele der durchgef ¨uhrten Arbeit sind:
• ein verbessertes Verst¨andnis von Mehrphasenstr¨omungen und der Prozesse in gekl ¨ufteten por¨osen Medien.
• die Entwicklung eines konzeptuellen Modells, das die Untersuchung einer Zweiphasenstr¨omung in Kl ¨uften beliebiger Gr¨oße, Ausrichtung und Form erlaubt.
• die Entwicklung eines flexiblen Kontinuumsansatzes, der die Vorteile des Kon-tinuumsansatzes mit denen diskreter Kluftmodelle kombiniert, sodass eine entsprechende ¨Ubertragung der Str¨omungseigenschaften individueller Kl ¨ufte auf effektive Parameter grobmaßst¨ablicher Bl¨ocke erm¨oglicht wird, w¨ahrend Str¨omungs- und Transporteigenschaften erhalten bleiben.
• die Entwicklung eines auf dem obigen konzeptuellen Modell basierenden effizienten und genauen numerischen Simulators.
• das Aufzeigen einer m¨oglichen Vorgehensweise zur Unterschung von Str¨omung in zerkl ¨ufteten por¨osen Medien.
Um die angestrebten Ziele dieser Arbeit zu erreichen, werden mehrere Schritte durchgef ¨uhrt:
1. Die Erarbeitung der Hauptgleichungen, welche Mehrphasenstr¨omungen in gekl ¨ufteten por¨osen Medien beschreiben, sowie das Verst¨andnis der Str¨omungsprozesse und der dominierenden Mechanismen. Auf einem fei-nen Maßstab, auf dem jede einzelne Kluft ber ¨ucksichtigt wird, wird eine ef-fiziente und zuverl¨assige Methode erstellt, die die Mehrphasenprozesse in Kluft-Matrix-Systemen anhand eines diskreten Kulftmodells (DFM) (Kap.2.5) beschreibt. Dieses Modell wird anschließend mittels des 2pDFM Moduls im Simulationsframework DuMuximplementiert (Kap. 3).
Im Gegensatz zu Einphasenmodellen spielen Kapillarkr¨afte in der Mehrpha-senstr¨omungen eine wichtige Rolle. Aus diesem Grund
a) werden grundlegende Beziehungen zwischen Kapillardruck und S¨attigung f ¨ur por¨ose Kluft-Matrix-Systeme hergeleitet, indem unterschied-liche Kluftverteilungen generiert werden und die S¨attigungsverteilung innerhalb der Systeme berechnet wird;
b) wird die Sensitivit¨at des nichtlinearen Verhaltens der Mehrphasen-str¨omungen in gekl ¨ufteten por¨osen Systemen bis zum Wechsel der grund-legenden Verh¨altnisse und der Grenzbedingungen mittels einer Sensiti-vit¨atsanalyse untersucht;
c) wird bei der Handhabung gekl ¨ufteter por¨oser Systeme eine ad¨aquate Be-schreibung der Grenzfl¨achen zwischen Kl ¨uften und Matrix notwendig. Dieses ist aufgrund der stark kontrastreichen hydraulischen Eigenschaf-ten der Kl ¨ufte und der Matrix, besonders im Hinblick auf die Mehrpha-senstr¨omung, eine Herausforderung. Aus diesem Grund m ¨ussen bereits vorhandene numerische Entw ¨urfe erweitert werden, um die Formulie-rungsschwierigkeiten der Grenzfl¨achen zu bew¨altigen.
2. Im n¨achsten Schritt geht es um die Entwicklung eines effektiven groß-maßst¨ablichen Kontinuumsmodells. Die Auswahl des Kontinuumsmodells erfolgt wegen des hohen Rechenaufwands diskreter Ans¨atze und wegen der in der Praxis fehlenden Informationen zu Details der Kluftgeometrien. Die Methode der geschachtelten interagierenden Kontinua (MINC) (in Kap. 2.7. im Detail beschrieben), die in die numerische Toolbox DuMuxim Modul 2pMINC
implementiert wird, stellt eine zuverl¨assige M¨oglichkeit dar. Die transiente Komponente wird mit ber ¨ucksichtigt durch die L¨osung der Dynamik innerhalb der Gesteinsmatrix.
3. Eine weitere Herausforderung ist die Bestimmung der effektiven hochskalierten Str¨omungsparameter (d.h. effektive Permeabilit¨at, effektive Por¨osit¨at, relative Permeabilit¨atsfunktionen, Mobilit¨atsgr¨oßen, usw.) und der Geometrien der Bl¨ocke eines grobmaßst¨ablichen Modells. Die MINC-Methode wird mittels eines Ansatzes erweitert, der sich der L¨osung einer Einphasenstr¨omung bei einem diskreten Feinskalenproblem bedient, um die Str¨omungsparameter eines grobmaßst¨ablichen Gitterblocks zu generieren.
4. Die vorgeschlagenen Methoden und die entsprechenden numerischen Simula-toren werden ¨uberpr ¨uft und getestet.
Allgemeine Zweiphasenstr ¨omungsgleichungen
Die Erhaltung der Masse im Kontext von Mehrphasenstr¨omung kann f ¨ur jedes Fluid in Phase α folgendermaßen ausgedr ¨uckt werden:
∂(φSαρα)
∂t + ∇ · (vαρα) − ραqα= 0 (1)
wobei die Geschwindigkeit vαmittels des erweiterten Darcygesetzes berechnet wird:
vα= −kµrα
αK · (∇pα− ραg) (2)
Daraus ergeben sich die allgemeinen Zweiphasenstr¨omungsgleichungen f ¨ur por¨ose Medien (Helmig [1997]):
∂(φSαρα)
∂t − ∇ · Kkrα
Gekl ¨uftete Systeme – ein ¨Uberblick
Entsprechend der Vielzahl der Skalen und der Schwierigkeiten, in denen Kl ¨ufte auftreten, werden die Kluftstr¨omungsmodelle grob in folgende drei Modelltypen ein-geteilt – diskrete Kluftmodelle (DFM), Kontinuumskluftmodelle (CFM) und hybride Modelle (s. Berkowitz [2002]).
DFM
Die DFM (so wie sie bei Bogdanov et al. [2003b,a], Hoteit and Firoozabadi [2005], Reichenberger et al. [2006], Assteerawatt [2008] genutzt werden) geh¨oren zur Nah-feldklasse (”near field”) der konzeptuellen Modelle und beschreiben die Effekte indi-vidueller Kl ¨ufte explizit. DFM stellen die genaueste Methode dar, um Str¨omungen in gekl ¨ufteten Formationen zu untersuchen. Ebenso k¨onnten sie auch eine Alternative zu den CFM darstellen, f ¨ur den Fall des Weglassens der Annahme, dass sich ein Kluft-system wie ein Kontinuum verh¨alt. Weiterhin k¨onnen hiermit die Effekte einzelner Kl ¨ufte explizit untersucht werden. Allerdings stoßen DFM bei der Anwendung auf Feldmaßstab auf zwei grunds¨atzliche Schwierigkeiten: sie ben¨otigen viel Rechenzeit und weisen sehr hohe Anforderung zur Generierung eines Kluftnetzwerks auf. Es gibt zwei M¨oglichkeiten Kl ¨ufte mittels DFM zu beschreiben:
• Der gleichdimensionale Ansatz (DFM-E): Kluft und Matrix werden mit Elemen-ten der gleichen Dimension diskretisiert. Dies bedeutet eine hohe Anforderung an die Netzwerkgenerierung sowie den numerischen Simulator (z.B. Gebauer et al. [2002]).
• Der niederdimensionale (gemischt-dimensionale) Ansatz (DFM-L) (z.B. Hel-mig [1997], Reichenberger et al. [2006]). Kl ¨ufte werden als niederdimensionale Elemente in Bezug auf die umgebende Matrix dargestellt; d.h. lineare (1D) Elemente in zweidimensionaler Umgebung und planare Elemente (2D) in drei-dimensionaler Umgebung.
DFM Annahmen
Die wichtigsten Annahmen der diskreter Kluftmodelle sind:
1. Die absolute Permeabilit¨at der Kl ¨ufte ist gr¨oßer als die absolute Permeabilit¨at der Gesteinsmatrix.
3. Die Mehrphasenstr¨omungsgleichungen (2.34) sind sowohl in der Gesteinsma-trix als auch im Kluftsystem g ¨ultig. Die Erweiterung des Darcygesetzes ist g ¨ultig und die Str¨omung ist f ¨ur alle Str¨omungsbereiche laminar.
4. Die relativen Permeabilit¨ats- und Kapillardruckfunktionen gelten sowohl f ¨ur Kl ¨ufte als auch f ¨ur Matrix. Der Kapillardruck wird als monoton abstei-gend angenommen und die Kapillardruckfunktionen f ¨ur Gestein und Matrix
¨uberkreuzen sich nicht.
5. Die benetzende Phase ist sowohl in Kl ¨uften als auch in der Gesteinsmatrix mobil.
Weitere notwendige Annahmen f ¨ur die Erweiterung der Schnittstellenbedingung f ¨ur Kluft-Matrix-Bereiche sind:
• Der Kapillardruck in der Matrix ist immer gr¨oßer als in den Kl ¨uften. F ¨ur die Brooks-Corey Formulierung bedeutet dies, dass der Eindringdruck in der Matrix gr¨oßer ist als derjenige in der Kluft.
• die Kontinuit¨at der Str¨omung an der Schnittstelle
• die Kontinuit¨at der intensiven Zustandsvariablen an der Schnittstelle: Kapillar-druckgleichgewichtsbedingung.
Die Gleichgewichtsgr¨oßen zwischen den Kluft- und den Matrixvariablen ergeben sich aus der erweiterten Kapillardruck-S¨attigung-Schnittstellenbedingung (mit der oberen Indizes M f ¨ur Matrix und F f ¨ur Kluft):
SM w = 1, if S F w> S∗w (pM c)−1 pFc SFw else (4) CFM
In großskaligen geologischen Formationen, g¨abe es Millionen Kl ¨ufte, die ber ¨ucksichtigt werden m ¨ussten. Dieses w¨are eine unglaubliche Aufgabe f ¨ur einen DFM Simulator. Daher wird in einem ¨ublichen CFM-Ansatz das Reservoir als zwei getrennte, interagierende Kontinua, n¨amlich Kluft und Matrix, aufgeteilt, die ¨uber einen Austauschterm (Transfer Term) gekoppelt sind. F ¨ur diesen Austauschterm wird Feindiskretisierung der DFM oder die detaillierte Kluftbeschreibung w¨ahrend der Simulation nicht ben¨otigt. Dieses sogenannte Doppelpor¨osit¨atsmodell (DP) wur-de von Barenblatt et al. [1960] and Warren and Root [1963] vorgeschlagen (Abb. 1 (2)).
Abbildung 1: 1) Gekl ¨uftetes Reservoir 2) schematische Repr¨asentation des MINC Modells (3D) 3) 2D Querschnitt durch das ”W ¨urfelzuckermodell”4) Konnektivit¨aten: a) Konnektivit¨at der grobmaßst¨ablichen Bl¨ocke oder Kluft-Kluft; b) Kluft-Matrix c) Matrix-Matrix.
Die Grenzen von CFM
Das klassische DP Modell st¨oßt schnell an seine Grenzen, wie z.B. in Pruess and Narasimhan [1982b], Zimmerman et al. [1996], Karimi-Fard et al. [2006], Unsal et al. [2009], Geiger et al. [2010] aufgezeigt.
Diese Methode funktioniert am besten, wenn das Kluftnetzwerk gen ¨ugend gut verbunden ist und der Unterschied zwischen den Permeabilit¨aten in Kl ¨uften und Matrix groß ist. Einige der Annahmen klassischer DP Modelle sind:
• die Str¨omung zwischen Kluft und Matrix, genannt Interpor¨osit¨atsstr¨omung, ist scheinkonstant (quasi-steady);
• die Str¨omung in der Matrix und die r¨aumliche Schwankung innerhalb der Matrixbl¨ocke wird vernachl¨assigt ;
• Verwendung empirischer Transferfunktionen zur Beschreibung der komplexen Kluft-Matrix Interaktionen;
• Verwendung von gemittelten Eigenschaften f ¨ur das Kluftnetzwerk, welche dessen Komplexit¨at in einem Permeabilit¨atstensor nicht erfassen k¨onnen;
• es existiert keine klare Prozedur zur Bestimmung der Parameter eines spe-zifischen diskreten Kluftsystems. Eine Reihe von Ans¨atzen versucht, verbes-serte Transferfunktionen zu bestimmen, die den Austausch zwischen Matrix und Kluft genauer darstellen [Zimmerman et al. [1993], Abushaikha and Gos-selin [2009]], oder die ¨aquivalenten effektiven Parameter zu quantifizieren (¨aquivalente Kluftpermeabilit¨at und Matrixblockdimensionen) [Zimmerman and Bodvarsson [1995],Min et al. [2004]].
Ein weiterer Ansatz zur Beschreibung der Matrixdynamik ist die Einf ¨uhrung eines Untergitters (subgridding) (z.B. Pruess and Narasimhan [1982b], Karimi-Fard et al. [2006]).
MINC
Die multiple interagierende Kontinua (MINC) Methode geh¨ort zur CFM Kategorie und ist eine Sch¨atzmethode zur Modellierung von Fluid- und Hitzestr¨omungen in gekl ¨ufteten por¨osen Medien. Entwickelt wurde sie von Pruess and Narasimhan [1982b]. Pruess [1992] enth¨alt einen kleinen Leitfaden zur MINC-Methode f ¨ur die Modellierung von Str¨omung und Transport in gekl ¨ufteten por¨osen Medien. Diese Me-thode ist eine Verallgemeinerung des Doppelpor¨osit¨atskonzeptes, jedoch behandelt sie die Interpor¨osit¨atsstr¨omung auf v¨ollig transiente Weise. Wenn das Kluftsystem gen ¨ugend verbunden ist, weisen einige der Bereiche kurz nach einer Perturbation fast die gleichen thermodynamischen Eigenschaften (Druck, S¨attigung, Temperatur, usw.) auf. Daher konnen diese Bereiche zu eigenen Kontinua zusammengefasst werden.¨ Der klassische MINC Ansatz nimmt an, dass die thermodynamischen Eigenschaften in der Matix je nach Distanz zur n¨achsten Kluft variieren. Mit anderen Worten – die ¨aquipotenziellen Fl¨achen haben den gleichen Abstand zur n¨achsten Kluft. Die-ses f ¨uhrt zur Aufteilung des Str¨omungsbereichs in berechenbare Volumenelemente, wobei alle Grenzfl¨achen/Schnittstellen der Volumenelemente parallel zur n¨achsten Kluft ausgerichtet sind, ¨ahnlich wie bei einer Matrioschkapuppe. Abb. 1 zeigt die Grundidee dieser Methode. Konzeptuell ¨ahnliche Modelle, bei denen die Str¨omung in den Matrixbl¨ocken auf einer Feinskala gel¨ost werden, wurden f ¨ur DP Modelle mittels asymptotischer Homogenisierung entwickelt (Douglas Jr. et al. [1991], Lewan-dowska et al. [2004]). Allerdings bleibt der Schwerpunkt der MINC Methode, wie bei jedem DP Modell, bei der Sch¨atzung der effektiven Permeabilit¨aten sowie, dar ¨uber hinaus, bei der Unterteilung des Matrixbereichs. Das klassische Konzept nutzt eine idealisierte Darstellung des gekl ¨ufteten Reservoirs (s. W ¨urfelzuckermodell). Die Interpor¨osit¨atsstr¨omungsparameter werden zu einem sp¨ateren Zeitpunkt kalibriert.
Weiterentwicklungen der MINC Methode
In den letzten Jahren ist die hydraulische und geologische Beschreibung gekl ¨ufteter Reservoire schnell vorangeschritten. Daher m ¨ussen die hochskalierten Parameter trozt der unbekannten Geometrie der Kl ¨ufte und auch f ¨ur den Fall, dass gen ¨ugend Informationen zur Kluftverteilung vorhanden sind, bestimmt werden. Wenn man diese vernachl¨assigt und stattdessen die gleichen Austauschterme (Transfer Term) f ¨ur das gesamte Gebiet nutzt, w¨aren die Resultate weniger genau. Das Ziel ist, sich die Genauigkeit der DFM und die rechentechnische Einfachheit der CFM zunutze zu machen. Karimi-Fard et al. [2006] f ¨uhren eine hochskalierende Methodik als Erweite-rung der klassischen MINC-Untergitterprozedur ein, die die effektiven Parameter von einer detaillierten Kluftbeschreibung ableitet (d.h. die Teilnetze werden f ¨ur jeden grobmaßst¨ablichen Block durch die Nutzung der Isodruck-Fl¨achen konstruiert, die man aus der lokalen Druckl¨osung einer DFM f ¨ur den ganzen Block erh¨alt). Einige der Hauptschwierigkeiten hinsichtlich der Hochskalierung gekl ¨ufteter Reservoire werden in Vitel and Souche [2007] zusammengefasst. Sie argumentieren, dass die auf die grobmaßst¨ablichen Bl¨ocke angewandten Grenzbedingungen das dynamische Verhalten in weiterer Entfernung des Blocks nicht ber ¨ucksichtigen.
Die Standard-Hochskalierungsans¨atze bilden Eindringdruckeffekte nur ungenau ab. Im Kontext eines sich neu ergebenden Anwendungsfeldes, der CO2Lagerung in
geologischen Formationen, erweitert diese Arbeit den von Karimi-Fard et al. [2006], Gong et al. [2008] untersuchten Anwendungsbereich, indem die Migration der nicht-benetzenden Phase in einen von der nicht-benetzenden Phase vollges¨attigten Bereich ber ¨ucksichtigt wird.
Die erweiterte MINC Methode
Die erweiterte MINC Methode ist eine Hochskalierungsmethode, die von Karimi-Fard et al. [2006] als eine Erweiterung der ehemaligen Teilnetzprozeduren entwickelt wurde (s. Pruess and Narasimhan [1982b]). Die Verbindungen zwischen dem Kluft-netzwerk und der Matrix sowie der des inneren eines grobmaßst¨ablichen Blocks werden anhand der L¨osung einer kompressiblen Einphasenstr¨omungsgleichung als DFM f ¨ur einen grobmaßst¨ablichen Block bestimmt:
∂(φρ) ∂t − ∇ ·
K
µ∇p = q. (5)
Die erweiterte MINC Methode beginnt mit einer detaillierten Beschreibung des Kluftnetzwerkes, welches mit einem DFM Ansatz zur Bestimmung der Bl¨ocke f ¨ur das Matrixstr¨omungsproblem gel¨ost wird. Dieses unterscheidet sich von fr ¨uheren Methoden, die eine idealisierte Darstellung der Matrix-Subregionen annehmen (s.
W ¨urfelzuckermodell, Abb. 1(2)) und daraus die Str¨omungsparameter kalibrierten. Die Methode besteht aus zwei Schritten: Im ersten Schritt werden die lokalen Teil-netze der grobmaßst¨ablichen Bl¨ocke bestimmt, welche die Matrixdynamik aufl¨osen. Durch die isolierte Betrachtung der grobmaßst¨ablichen Bl¨ocke wird der lokale Fluid-austausch zwischen dem Kluftnetzwerk und der Matrix sowie der geschachtelten Matrixvolumina untereinander berechnet. Gleichzeitig werden die geometrischen Parameter (Großporenvolumina) berechnet (Abb. 2).
Im zweiten Schritt wird die hochskalierte Zwischenblockstransmissibilit¨at mittels einer Zweipunkt-Sch¨atzmethode bestimmt, indem eine permeameter¨ahnliche Grenz-bedingung f ¨ur jeden grobmaßst¨ablichen Nachbarblock (Ωi, Ωj) festgesetzt wird und
erneut eine stabile Einphasenstr¨omungsgleichung gel¨ost wird (Abb. 3). Ti, j= Qi, jµ
ρ pi− pj (6)
Das mathematische Modell f ¨ur die MINC Methode
Das grobmaßst¨abliche Modell wird durch Gleichungen allgemeiner Form beschrie-ben, wie der Massenerhaltungsgleichung f ¨ur Kluftkontinua (F):
∂(φSαρα)F
∂t − ∇FF· Ke f fk
rα
µα ρα(∇pα− ραg) F− ραqα− τ = 0, (7)
und f ¨ur Matrixkontinua Mjmit j = 1,J:
∂(φSαρα)Mj ∂t − ∇MM· Ke f fk rα µα ρα(∇pα− ραg) Mj = 0, (8)
Abbildung 2: Darstellung des Druckprofils und der Isodrucklinien in quasi-station¨are Zustand f ¨ur den lokalen Block Ωiauf der Feinskala und Konstruktion
der idealisierten verschachtelten Volumenelemente Ωi
1,Ωi2,...,Ωi6der
Abbildung 3: Bestimmung der grobmaßst¨ablichen Zwischenblocktransmissibilit¨at der grobmaßst¨ablicher Bl¨ocke Ωiund Ωjmit mittlerem Blockdruck pi
und pj τ = −1 V Vm ∂(φSαρα) ∂t d(MM), (9)
wobei ∇FFden Grobskalenoperator (Wechselwirkungen zwischen den Kluftbl¨ocken,
die auf globaler Skala gel¨ost werden), und ∇MMden Feinskalenoperator
(Wechsel-wirkungen innerhalb der Grobskalenbl¨ocke zwischen den Matrixkontinua, die auf einer lokalen Skala gel¨ost werden) anzeigt.
Der Austauschterm (Transfer Term) τ modeliert den Austausch zwischen der Matrix und der Kluft und wird von der L¨osung der lokalen Matrixstr¨omungsgleichung bestimmt.
(Ke f f)
Mjstellt die effektive (hochskalierte) Matrixpermeabilit¨at der Kontinua Mjdar.
Numerische Implementierung
Numerische Modellierung
Die Erhaltungsgleichungen, die im mathematischen Modell beschrieben wurden, bilden ein System stark gekoppelter, nichtlinearer partieller Differentialgleichungen. Diese Systeme sind zu komplex, um analytisch gel¨ost zu werden und werden deshalb numerisch gel¨ost. Das System partieller Differentialgleichungen wird r¨aumlich mit der BOX-Methode und zeitlich mit dem vollimpliziten Eulerschema diskretisiert. Die Linearisierung wird mit der Newton Methode durchgef ¨uhrt. Danach werden die Zweiphasenkluftstr¨omungsmodelle, das DFM und das MINC Modell im nu-merischen Framework DuMuxdurch die Module 2pDFM und 2pMINC implementiert.
Damit geh¨oren sie zur sogenannten Kategorie erweiterter Modelle (Flemisch et al. [2011a]).
Ablauf (Workflow)
Eine Simulation kann in drei Abschnitte unterteilt werden: Pre-Processing, Proces-sing und Post-ProcesProces-sing. Im Pre-ProcesProces-sing gibt es zwei Hauptstufen. Die erste ist die Generierung des diskreten strukturellen Modells des Kluftbereichs, die Konver-tierung der Geometriedaten mittels eines Gitternetzgenerator und die Erstellung des Gitters. Im zweiten Schritt werden die Kl ¨ufte mittels des Programms artreader berechnet sowie die Grobmaßstabsblockparameter bestimmt.
Der Processing-Abschnitt berechnet die L¨osung eines gegebenen Str¨omungsproblems. Die Module 2p, 2pDFM und 2pMINC basieren hierbei auf einem vollimpliziten, voll-gekoppelten Diskretisierungsverfahren. Die geometrische Modellierung und die Gitternetzerstellung im Pre-Processing, das Processing sowie die Visualisierung der Ergebnisse werden lose an die DuMuxSimulationsumgebung gekoppelt.
Modellvalidierung
Die numerischen Modelle werden validiert, indem numerische Resultate mit experi-mentellen Daten sowie weiteren zuverl¨assigen Simulatoren verglichen werden.
DFM: Es wurden mehrere Experimente durchgef ¨uhrt, um das DFM bez ¨uglich der
Fluidtypen, der nieder- und gleich-dimensionalen Darstellung der Kl ¨ufte und der Ausrichtung der Kl ¨ufte zu testen. Der DFM Zweiphasensimulator wird mit anderen Simulatoren verglichen, s. Wu et al. [2004b] , wie auch mit den Laborexperimenten von Kazemi [1976]. Kazemi f ¨uhrte eine Reihe von Laborexperimente zur Wasserim-bibition gekl ¨ufteter Matrixkerne unter Verdr¨angung von ¨Ol durch. Ein großer Teil der Modellvalidierung bestand darin, aufzuzeigen, dass die Simulationsergebnisse der niederdimensionalen DFM mit denen zuverl¨assiger gleichdimensionaler DFM ¨ubereinstimmen. Daher werden diese als Referenzl¨osung zur Validierung des MINC Modells eingesetzt. Ein weiteres untersuchtes Ph¨anomen ist der Einfluss der Feinheit des Gitters auf die L¨osung.
MINC: F ¨ur das MINC Modell wurden sechs Validierungsbeispiele erstellt, die
be-weisen, dass das Modell zuverl¨assig ist.
Anwendung des Modells auf der Feldskala
Das erweiterte MINC Modell wird auf mehreren Feldskalenprobleme angewendet. Der erste Test ist ein idealisiertes periodisches Kluftsystem. Die effektiven Para-meter der erweiterten MINC Methode sowie die Referenzl¨osung werden mit einer
DFM berechnet. Die Tests untersuchen den Einfluss der Str¨omungsrichtung und der Kluftorientierung auf die S¨attigungsverteilungs- und die Frontausbreitungsge-schwindigkeit. Im ersten Test, einem permeameter¨ahnlichen Experiment, verl¨auft die Str¨omung parallel zur Ausrichtung der horizontalen Kl ¨ufte, wobei diese im zweiten Experiment, einem ”Quarter Five SpotProblem, diagonal zur Kluftausrich-tung verlaufen. Beide Beispiele eine gute N¨aherung der S¨attigungsverteilungs- und Frontausbreitungsgeschwindigkeit. Die MINC Methode erlaubt die Vorhersage der S¨attigung innerhalb des Matrixbl¨ocke, wenn der Eindringsdruck ¨uberwunden ist. Das n¨achste Beispiel ist ein nat ¨urlich gekl ¨uftetes Reservoir in Bristol mit einem stark verbundenen Kluftsystem. Die Tests zeigen, dass der 2pMINC Simulator die Str¨omungen und S¨attigungsprofile, die sich aus dem diskreten 2pDFM Simulator er-geben, gut approximieren kann und auch sehr hohe Rechenzeitersparnisse ergibt. Gleichzeitig werden sowohl der Matrix-Kluft-Transfer wie auch die Eindringsdruck-effekte korrekt gesch¨atzt; sie gr ¨unden auf einem physikalisch korrekten Ansatz. Die Resultate werden als Profile der nicht-benetzenden Phase entlang der Diagonale (von links unten nach rechts oben) dargestellt (Abb. 4).
Abbildung 4: S¨attigungsprofile der nicht-benetzenden Phase auf der Diagonalen f ¨ur das permeameter¨ahnliche Experiment. Obere Werte beziehen sich auf DFM L¨osungen, untere Werte wurden mittels des MINC Modells berechnet.
Abschließende Anmerkungen
• Aus dem niederdimensionalen DFM ergeben sich Systeme, die viel einfacher gel¨ost werden k¨onnen als jene der gleichdimensionalen DFM. Das DFM-L ben¨otigt kein Teilnetz innerhalb der Kluft und vermeidet daher kleine Elemente innerhalb der Kl ¨ufte.
• Die BOX Methode beachtet das Gesetz der Massenerhaltung und kann pro-blemlos auf unstrukturierte Gitternetze angewendet werden.
• Die MINC Methode kann den Matrixeinfluss mitber ¨ucksichtigen, der eine große Rolle bei großskaligen Mehrphasenstr¨omungsprozessen spielt. Dieses ist wichtig f ¨ur die Evaluation der Speicherkapazit¨at und des Reservoirverhaltens. • Die erweiterte MINC Methode ist flexibel, da sie sowohl die Genauigkeit der DFM als auch die niedrige Rechenzeit der CFM vereint. Dieses ist ein Kom-promiss zwischen der Genauigkeit und der Rechengeschwindigkeit, der dem Modellierer ¨uberlassen wird. Die Kl ¨ufte werden auf rechnerisch effizienten Gitternetzen abgebildet.
• Die erweiterte MINC Simulationsmethode ist allgemein anwendbar, da sie die Bew¨altigung sowohl dicht gekl ¨ufteter Systeme mit unregelm¨aßigen Geo-mentrien erlaubt, als auch von Topologien, deren Kl ¨ufte schw¨acher verbunden sind.
• In der erweiterten MINC Methode bleiben die gesamten Str¨omungs- und Trans-porteigenschaften erhalten.
1.1 Motivation
A fracture is considered a separation of an object or material as a result of mechanical failure. Fractures may occur in every natural or synthetic material, in biological tissues like bones and even in the best engineered materials. For instance, fractures can be encountered in concrete blocks, bricks used as structural material for con-structions, geosynthetic clay liners used in landfills and disposal pits, underground storage tanks, asphalt of the roads, ceramic plates, in enamel of the tooth, geologic formations etc (Figure 1.1 and 1.2).
Fractured systems play an important role in the fluid flow and transport in porous media, while many engineering applications require a correct understanding of these processes. The main focus of this research is on the fractured geological formations. As indicated by Berkowitz [2002] “fractured geological formations are ubiquitous throughout the world, and are of interest in a number of contexts :
1. reservoir exploitation for water supply;
2. contamination from subsurface waste repositories; 3. petroleum reservoir exploitation;
4. geothermal reservoir exploitation and heat storage;
5. mining and mineralization processes (in situ leaching and location of ore bod-ies);
6. geotechnical applications (including effects on underground storage reservoirs, tunnels and other structures); and
7. deeper Earth systems, such as earthquakes and ocean floor hydrothermal venting.”
An emerging field of application addresses to the mitigation of the anthropogenic effects on climate change through geological sequestration of carbon dioxide (IPCC [2005], Klara et al. [2003], Bachu [2003], Kopp et al. [2009]). Other line of research is investigating possibilities to dispose the radioactive waste in rock dumping sites,
(a) Dry cracked earth (www.123rf.com) (b) Heavily fractured exposed rock (Keller [1996])
(c) Fractured asphalt (www.123rf.com) (d) Fractured asphalt
(e) Fractured limestone with faults from the
Bristol Channel, UK (csmp.ese.ic.ac.uk) (f) Fracture corridor in a clastic rock, Algeria(Iain Bush, Schlumberger)
(g) Electron microscope scan of a bone
(Mas-trolorenzo et al. [2010]) (h) Scanning electron microscopy micrographof an ELAT-DS electrode in a fuel cell (Acosta et al. [2006]).
(a) Microfissures in the Japanese ceramics
-Raku (ceweb2.uml.edu) (b) Japanese ceramics (Raku): The crackle pat-tern is caused by carbon being trapped in tiny fissures that heal over as the glaze cools (www.123rf.com)
(c) Skull fracture (www.umm.edu) (d) Fractured-like structure of human brain (www.joelertola.com)
(e) Fractured icy surface on Europa, Jupiter’s
satellite (explanet.info) (f) Cemented fractures in mountain insideGale crater on Mars, taken by NASA’s Mars Reconnaissance Orbiter (mars.jpl.nasa.gov)
as they are considered virtually impermeable (e.g. Bodvarsson et al. [2003], Re-ichenberger et al. [2003]). For this kind of applications are of interest the degassing effects that may result due to pressure drops in fractures in the vicinity of deep high-radioactive-waste-disposal sites.
Ideally, the barriers of the storage sites have to be impermeable so that the contained substances can not escape. However, because these barriers are usually interlaced with fractures and shear zones which are generally interconnected, the clear under-standing of multiphase flow processes in porous and fractured porous media is of capital importance. In fact, “all the geological media present a certain degree of fracturing” (Neuman [2005]).
Computer modeling is the process of developing mathematical representations of a real-life or hypothetical situations so that they can be better studied and to gain insight on how those systems function and could be operated. The subsurface waste-disposal sites have to guarantee the safety for hundreds, if not thousands of years, therefore the key to understanding the multiphase processes and foresee the behavior of the reservoir is to construct an accurate model. Some of the greatest challenges for modeling multiphase flow in fractured porous media are caused by the big discrepancies in the characteristics of the fractures and the surrounding rocks and the difference in the time-scales in flow processes.
In this research mathematical models and numerical simulation tools are developed for fractured porous media that can play an important role to evaluate the feasibility of the disposal subsurface reservoirs, to designing and analyzing field tests and to improve and optimize the operation of those geologic disposal systems.
1.2 Objectives of the Research
The main objectives of this study are to:
• improve the understanding over the concepts of the multiphase flow and processes in the fractured porous media
• develop a conceptual model that allows the study of two-phase flow in fractures of arbitrary size, orientation and shape.
• develop a flexible continuum approach that combines the advantages of the con-tinuum and discrete fracture model concepts so that it allows to appropriately transfer the flow characteristics of individual fractures to effective parameters assigned to coarse blocks while preserving the flow and transport features. The flexibility in this context refers to the ability of choosing the accuracy of the solution, the computation speed, and the ability of working with spatial information about the fracture system of various complexity and detail.
• develop an efficient, accurate numerical simulator based on the above men-tioned conceptual model
• construct a workflow for studying flow in fractured porous media.
1.3 Thesis Outline
Several steps have to be implemented for achieving the imposed objectives of the research:
1. Establishing the governing equations that describe the multiphase fluid flow in fractured porous systems and understanding the flow processes together with the dominating mechanisms. At the fine scale, which considers each individual fracture, an efficient and reliable method is set up for describing the multiphase processes in fracture-matrix systems through a discrete fracture model (DFM) (Sec. 2.5). The solution is implemented in the numerical simulator 2pDFM (Sec. 3) under the general frame program DuMux.
In contrast to single-phase models the capillary forces play an important role when coming to multiphase flow, therefore:
a) constitutive relationships between capillary pressure and saturation for fracture-matrix porous systems are determined by generating different fracture aperture distributions and calculating the saturation distribution within the systems (Sec. 4.1.2) ;
b) the sensitivity of the non-linear behavior of the multiphase flow in frac-tured porous system to the change of constitutive relationships and bound-ary conditions is investigated by means of sensitivity analysis;
c) dealing with fractured porous systems requires an adequate description of the interface between fractures and matrix. This is a challenging task, especially concerning multiphase flow, due to the strong contrasts in hydraulic properties of fractures and matrix. For this reason existing numerical schemes must be extended to be able to cope with interface formulation difficulties.
2. The next step is to develop an efficient large-scale continuum model (Figure 2.6). The choice of a continuum model is motivated by the high computational demands of a discrete approach and by the practical reason of lacking detailed information about the fracture geometries when dealing with real field applica-tions. The multiple interacting continua (MINC) method (described in detail in Sec. 2.7), implemented in the numerical toolbox DuMuxunder the name
2pMINC simulator (Ch.3), represents a trustworthy option because it takes the transient component into consideration by resolving the dynamics within the rock matrix.
3. A challenging step is the determination of the effective upscaled flow pa-rameters (i.e. effective permeability, effective porosity, relative permeability functions, mobility terms, etc.), and the coarse-block geometries of the coarse model. The MINC method is extended with an approach that uses the solution of a single phase flow on a discrete fine-scale problem to generate the flow parameters for the coarse-scale grid block.
4. The proposed methods and the corresponding numerical simulators have been verified1, validated2and tested.
1.4 Structure
The thesis has the following structure:
• Chapter 2 deals with the development of the conceptual model. The equations that govern the processes in fractured porous media are introduced and an overview over the discrete and the continuum fracture model concepts (Sec. 2.5, Sec. 2.6) is given.
• Chapter 3 discusses technical aspects of the implementation of MINC in the 2pMINC simulator and DFM into 2pDFM. Also, it presents the developed work-flow for modeling two-phase work-flow in fractured porous media.
• In Chapter 4, the two fracture flow models, DFM and MINC, are verified, validated and tested by comparing simulation results with the experiments in the literature or to trusted reference solutions. The DFM examples enumerate: the comparison to laboratory results, the comparison of the equi-to lower-dimensional fracture approaches. The MINC verification of the inter-block connectivity is done with: an analogous single homogeneous block, a single homogeneous block with a heterogeneous interior and a quasi one-dimensional domain; while the coarse(-MINC)-block- and inner -connectivity is verified with an analogy to a two-dimensional domain in a series of nested blocks. Further investigations included the influence of the entry pressure for the viscous and pressure driven flows.
1In computer modeling and simulation, verification is the process of determining that a model
or simulation implementation accurately represents the developer’s conceptual description and specifications.
2validation is the process of determining the degree to which a model or simulation is an
accu-rate representation of the real world from the perspective of the intended uses of the model or simulation.
• Chapter 5 introduces numerical simulation examples performed with both ap-proaches, DFM and MINC, for an idealized fracture system and a real naturally fractured reservoir. The simulation results are compared in terms of accuracy and computational time.
• Chapter 6 summarizes the thesis with conclusions and possible future develop-ments.
Media
This chapter introduces the theory that stays behind the developed fracture models. It starts with an introduction into the fundamental definitions and concepts of porous media, then it gives a short classification of the fracture flow models and after that it presents their attributes and the limitations that lead to further improvements. Later are introduced the governing equations of the two-phase flow and the numerical models.
2.1 Scales of Consideration
Molecular scale [∼ 10−10m] is the scale where the individual movement of the
molecules and their interactions can be described. In chemistry and physics, the Avogadro constant indicates that each mole of substance has about 6.02214 × 1023
molecules. The computational cost to model just a mole of substance would be tremendous.
Microscale [∼ 10−6− 10−2m] is the scale derived from molecular scale averaging.
At this scale the matter appears to be continuous and the molecular movement is blended to the observer. The interfaces between fluids and/or solids are clearly defined. Also the microscale implies that the discontinuities (e.g. micro fissures, cracks), which are small compared to the pore size diameter or the fracture aperture, can still be recognized. New state variables are defined (e.g. viscosity, density), which are continuous within the space and can be differentiable functions of space and time. This consideration of continuum is a fundamental concept of fluid mechanics which allows the derivation of the balance equations (e.g. Navier-Stokes).
Macroscale [∼ 10−3− 10−1m] is the scale derived from averaging the microscale
properties within a representative elementary volume (REV). The averaging volume needs to be large enough to ensure that the average properties are independent of the size of the REV (Figure 2.1) but it should be much smaller than the size of the entire domain cf. Bear [1988]. The transition from the microscale to macroscale leads to new basic equations (e.g. Darcy’s law (see Sec. 2.3.1) ) with new parameters (e.g. saturation, permeability).
