Dissertation_AndreasKopp
Dissertation_AndreasKopp
Dissertation_AndreasKopp
Vorgelegt von
Andreas Kopp
aus Furtwangen
2009
Heft 182 Evaluation of CO2 Injection
Processes in Geological
Formations for Site Screening
von
Dr.-Ing.
Andreas Kopp
Kopp, Andreas:
Evaluation of CO2 Injection Processes in Geological Formations for Site Screening – / von Andreas
Kopp. Institut für Wasserbau, Universität Stuttgart. - Stuttgart: Inst. für Wasserbau der Univ.,
2009
Herausgegeben 2009 vom Eigenverlag des Instituts für Wasserbau, Universität Stuttgart
Druck: Document Center S. Kästl, Ostfildern
Danksagung
Ich möchte mich an dieser Stelle bei all den Personen bedanken, die zum Gelingen dieser
Arbeit beigetagen haben.
Zuerst möchte ich mich bei meinem Hauptberichter Rainer Helmig bedanken. Er gab mir
die Chance nach einjähriger beruflicher Tätigkeit an die Universität zurückzukehren und
die Vorzüge einer unabhängigen Forschungsarbeit zu genießen. Die große Freiheit die er mir
dabei gewährte und die konsequente Rückendeckung die ich permanent erfahren durfte waren
ein großer Ansporn für mich und Grundlage für die Freude bei der Arbeit. Seine fundierte
Fachkenntnis, die zahlreichen Kontakte zu anderen Forschungseinrichtungen und Institutio-
nen die er mir vermittelte und die fortwährende Leidenschaft und Motivation die fachlichen
Probleme zu lösen, haben diese Arbeit maßgeblich positiv beeinflusst. Bei Prof. Christoph
Clauser (Aachen) und Prof. Helge Dahle (Bergen) möchte ich mich für die Übernahme der
Mitberichte bedanken.
Mein besonderer Dank gilt Holger Class, der mir ständig mit fachlichen Ratschlägen und
Problemlösungen den richtigen Weg gezeigt hat. Seine Ruhe und Besonnenheit in schwieri-
gen Zeiten waren Grundpfeiler für die Durchführung der Arbeit. Unsere zahlreichen Reisen
zu Konferenzen, Workshops, Tagungen, etc. waren nicht nur lehrreich, sondern haben auch
viel Spaß gemacht.
Herzlicher Dank geht an meine Mitarbeiter am Institut für Wasserbau und insbesondere an
die Kollegen am Lehrstuhl für Hydromechanik und Hydrosystemmodellierung. Es herrschte
immer eine sehr gute Stimmung die sich teilweise im Privaten fortsetzte.
Zu guter Letzt möchte ich mich bei meiner Familie bedanken, die mich stets unterstützt hat
und ohne die ich diese Ausbildung nicht geschafft hätte. Meinem Onkel Siegfried möchte ich
besonders danken für stete Unterstützung und die finanzielle Hilfe während meines Studiums.
Meiner Lebensgefährtin Andrea möchte ich für den Rückhalt und das Verständnis in den
letzten Jahren und für unsere wundervolle Tochter Milena danken.
All religions, arts and sciences are branches of the same tree.
All these aspirations are directed toward ennobling man’s life,
lifting it from the sphere of mere physical existence
and leading the individual towards freedom.
Albert Einstein
German-Swiss-American physicist
1879 - 1955
Contents
Abstract I
Zusammenfassung i
1 Introduction 1
1.1 Trapping Mechanisms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.2 Objective of the Study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.3 State of the Art . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
I
II Contents
4 Dimensional Analysis 51
4.1 Derivation of Dimensionless Formulation . . . . . . . . . . . . . . . . . . . . 51
4.1.1 Fractional Flow Formulation . . . . . . . . . . . . . . . . . . . . . . . 51
4.1.2 Characteristic Values . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
4.1.3 Dimensionless Numbers . . . . . . . . . . . . . . . . . . . . . . . . . 54
4.1.4 Dimensionless Pressure and Saturation Equations . . . . . . . . . . . 55
4.2 Preliminary Definitions of Risk and Storage Capacity . . . . . . . . . . . . . 56
4.3 Analytical Investigations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
4.3.1 Dimensionless Numbers Dependent on Selections of Characteristic Values 57
4.3.2 Dimensionless Functions A,B, and C . . . . . . . . . . . . . . . . . . 60
4.3.3 Dimensionless Gradients . . . . . . . . . . . . . . . . . . . . . . . . . 61
4.4 Numerical Investigations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
4.4.1 Plume evolution in a 1-D gravitation-free reservoir . . . . . . . . . . . 62
4.4.2 Plume evolution in a radially symmetric 3-D reservoir . . . . . . . . . 66
4.4.3 Qualitative Dependencies of Risk and Storage Capacity on Dimension-
less Numbers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
4.5 Summary and Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
6 Sensitivity Analysis 87
6.1 Discussion of Sensitivity Analysis . . . . . . . . . . . . . . . . . . . . . . . . 87
6.1.1 The Morris Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
6.2 Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
6.3 Numerical Investigations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
Contents III
7 Risk Analysis 99
7.1 Discussion of Risk . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
7.1.1 Risk Scenarios in CO2 storage . . . . . . . . . . . . . . . . . . . . . . 100
7.1.2 Time aspect . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
7.1.3 Screening and Ranking Frameworks . . . . . . . . . . . . . . . . . . . 101
7.1.4 Risk Analysis Concept . . . . . . . . . . . . . . . . . . . . . . . . . . 101
7.2 Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
7.2.1 Primary Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
7.2.2 Secondary Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . 107
7.2.3 Procedure of Defining a Simulation Case . . . . . . . . . . . . . . . . 111
7.3 Numerical Investigations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112
7.3.1 Model Set-up . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112
7.3.2 The “CO2 Community Grid” . . . . . . . . . . . . . . . . . . . . . . 112
7.3.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
7.4 Qualitative Sensitivity Considerations . . . . . . . . . . . . . . . . . . . . . . 119
7.5 Summary and Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122
Bibliography 133
C Tables 161
3.1 Histograms data of reservoir parameters derived from the NPC database. . . 44
3.2 Kolmogorov-Smirnov test on absolute permeability and geothermal gradient. 46
3.3 Test on mutual reservoir parameter interrelations. . . . . . . . . . . . . . . . 47
4.1 Variation of Gravitational Number and Capillary Number versus the charac-
teristic velocity. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
4.2 Gravitational Number versus Capillary Number for a varying characteristic
velocity. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
4.3 Functions A, B, and C for a Brooks & Corey relative permeability relation
model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
4.4 Carbon dioxide saturations in a horizontal 1-D reservoir after 4 years modeltime. 62
V
VI List of Figures
7.1 Sketch of the radially symmetric model domain and definition of the leakage. 103
7.2 Histogram of relative frequency of permeability anisotropy derived from a model.108
7.3 Correlation functions between porosity and absolute permeability and NPC
database values. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
7.4 Capillary pressure dependence on water-rich phase saturation and porosity. . 111
7.5 Slice of the radially symmetric domain showing CO2 -rich phase saturation for
two cases. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
7.6 Calculated likelihood of failure surface. . . . . . . . . . . . . . . . . . . . . . 114
7.7 Calculated consequence of failure surface. . . . . . . . . . . . . . . . . . . . . 115
7.8 Calculated risk surface. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116
7.9 Calculated risk versus time for selected leaky well distances. . . . . . . . . . 117
7.10 Comparison with a literature leakage rate versus time. . . . . . . . . . . . . 121
A.1 Basis function (N) for the respective node in the 1-D case. . . . . . . . . . . 153
A.2 Finite element and finite volume mesh. . . . . . . . . . . . . . . . . . . . . . 155
A.3 Weighting function for the respective node in the 1-D case. . . . . . . . . . . 156
List of Tables
6.1 Parameters investigated in the sensitivity analysis, parameter range, and lit-
erature source. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
6.2 Qualitative ranking of parameter effect as the average of individual overall
parameter effects. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
A.1 Primary variables and substitution criteria for the 2p2cni-module. . . . . . . 152
C.1 Power-fitted coefficients A and B in Equation 7.7 to calculate risk contour lines.161
C.2 Coefficients Ai fitted to Equation 7.8 to calculate time contour lines. . . . . . 162
VII
Nomenclature
The following table shows the symbols used in this thesis. Local notations are explained in
the text. Vectors and matrices are written in bold type.
IX
X Nomenclature
subscripts
CO2 CO2 -rich phase
cr characteristic value
crit critical, at critical point
heat referring to heat
mass referring to mass
n non-wetting phase
pm porous medium
s solid phase
w water-rich phase
α, β referring to phases α, β
superscripts
C component
CO2 CO2 component
NaCl salt component
r trajectory
t at tangential point of fractional flow function
w water component
Abstract
The concentration of greenhouse gases in the atmosphere has increased due to tremendous
human fossil fuel consumption since the Industrial Revolution. This is most likely the cause
for an observed global increase in the average temperature and for the changing climate. It
is expected that further global warming will have drastic ecological and economic impacts.
No single technology will be sufficient to achieve the necessary emission reductions. Carbon
dioxide capture and storage (CCS) is a promising technology which could make a substantial
contribution. It is a process which captures CO2 from large local sources and then stores
it away from the atmosphere. Storage capacity estimates for deep saline aquifers are most
promising. The initial procedure for selecting a few aquifers for a CCS project is called site
screening. Aquifers identified in site screening then have to prove their suitability in further
investigations. Site screening is a challenging task, since usually few data are available and
the prognosis of the complex processes occurring in a reservoir after CO2 injection is difficult.
This study aims at improving the insight into CO2 injection processes in geological forma-
tions to assist site screening. The criteria in site screening include the estimation of the
storage capacity, which should be sufficient to store the long-term production of the CO2
source, and the long-term ability to store CO2 , which is related to the efficiency of the project
and risk arising due to possible CO2 leakages.
At first, the statistical characteristics of storage sites in potential geological formations are
calculated by analysis of a large database. The parameter ranges and distributions are used to
define typical reservoirs and serve as a basis for generating random reservoir setups respecting
statistical characteristics. The relation of forces occurring in reservoirs after CO2 injection
is analysed by dimensional analysis. By the identification of dominant forces and processes,
reservoirs with different parameter setups are compared with respect to their potential CO2
storage capacity and risk. A sophisticated concept for estimating the CO2 storage capacity
of geological formations is developed. Detailed, time-dependent storage-capacity estimates
are calculated in numerical experiments. The results are interpreted using the simultane-
ously calculated ratios of forces. The influence of individual reservoir parameters on storage
capacity and risk is investigated in a sensitivity analysis. Finally, a risk analysis on potential
CO2 leakage through pre-existing wells is performed. In numerous numerical experiments,
individual parameters are randomly sampled from the statistical parameter distributions
and leakage is calculated. A risk surface is derived which represents the average risk for CO2
leakage through pre-existing wells for any site with unknown reservoir properties.
I
Zusammenfassung
Aufgrund des enormen Verbrauchs an fossilen Brennstoffen in den letzten 160 Jahren,
stieg die Konzentration der Treibhausgase in der Atmosphäre stark an. Dieser Anstieg
der Treibhausgaskonzentrationen ist mit größter Wahrscheinlichkeit die Ursache für den
weltweiten Temperaturanstieg und die beobachteten Klimaveränderungen. Man erwartet,
dass ein weiterer Temperaturanstieg zu tiefgreifenden ökologischen Veränderungen und
ökonomischen Belastungen führen wird. Eine einzelne Technologie oder Maßnahme wird
die nötige Verringerung der Treibhausgasemissionen nicht leisten können, deshalb muss
gleichzeitig eine ganze Reihe an Maßnahmen ergriffen werden. Zu diesen Maßnahmen
gehören z.B. eine effizientere Energiegewinnung und -nutzung, der Ausbau der Nutzung
regenerativer Energien, die erhöhte Verwendung treibhausgasarmer Brennstoffe, die Ab-
scheidung des CO2 im Abgasstrom von großen CO2 -Produzenten und die anschließende
Einlagerung in tiefe geologische Schichten oder der Tiefsee (CCS), eine Umstellung auf eine
weniger klimabelastende Landwirtschaft, die Aufforstung, sowie der eventuelle Ausbau der
Kernenergie zur Stromerzeugung.
Der Fokus dieser Arbeit liegt auf der Abscheidung und Speicherung von CO2 in tiefen
geologischen Schichten. Das CO2 wird hierbei im Abgasstrom von z.B. Kohle, Gas
oder Öl betriebenen Kraftwerken abgeschieden, mittels Rohrleitungen zur Speicherstätte
geleitet und dann in eine mindestens 800 m unter der Tagesoberfläche liegende, poröse,
mit Salzwasser gefüllte Gesteinsschicht eingepresst. Nach der Injektion breitet sich das
CO2 als freie Phase lateral aus. Gleichzeitig steigt es aufgrund der wesentlich geringeren
Dichte gegenüber dem Salzwasser auf. Ein gewisser Anteil löst sich dabei im Salzwasser.
Um ein weiteres Aufsteigen zu verhindern, werden Formationen für ein solches Vorhaben
ausgewählt, die über eine geringdurchlässige Deckschicht (Caprock) verfügen. An dieser
Deckschicht wird das CO2 durch den nicht zu überwindenden Eindringdruck aufgehalten.
Die Deckschicht kann jedoch auch geologische Schwachstellen aufweisen oder kann durch
menschliche Aktivitäten, z.B. gebohrte Brunnen, potentielle Fließpfade für ein Entweichen
des CO2 aufweisen. Dies verursacht zum einen das Risiko für ein Scheitern des Vorhabens,
nämlich das Abtrennen des CO2 für geologische Zeiträume aus dem Kohlenstoffkreislauf
und der Atmosphäre, zum anderen besteht ein Gesundheitsrisiko von Menschen die mit
ausgetretenem CO2 in Kontakt kommen sowie der Beeinträchtigung von Tieren, Pflanzen
und ganzer Ökosysteme.
i
ii Zusammenfassung
Der Projektzyklus eines typischen CO2 Speichervorhabens beinhaltet zunächst die Ab-
schätzung der zu erwartenden Emissionen über den gesamten Produktionszeitraum des
Kraftwerks (o.ä.). Anschließend erfolgt die Suche nach geeigneten Formationen in der
Region. Wurde eine Vorauswahl getroffen, erfolgt für die identifizierten Standorte eine
detaillierte Charakterisierung durch die Entwicklung eines geologischen Modells, die Unter-
suchung der CO2 Ausbreitung und des Druckaufbaus durch numerische Simulationen, sowie
eine umfassende Risikobewertung. Letztendlich wird ein geeigneter und sicherer Standort
ausgewählt und ein Antrag über eine CO2 Speicherung an die zuständige Genehmigungs-
behörde gestellt. Nach erfolgter Genehmigung wird dann für bis zu mehreren Jahrzehnten
CO2 injiziert sowie die Ausbreitung des CO2 und der Druckaufbau innerhalb der Formation
überwacht. Abschließend wird der Standort stillgelegt und weiter überwacht.
Mit der vorliegenden Arbeit soll das Prozessverständnis von CO2 Injektionen in geologische
Formationen verbessert werden um die anfängliche Standortauswahl innerhalb einer Region
zu unterstützen. Diese Phase eines Projekts ist typischerweise durch einen Mangel an
detaillierten Standortinformationen gekennzeichnet. Um eine Vorauswahl treffen zu können,
muss dennoch die Speicherkapazität einer Formation abgeschätzt werden. Außerdem sollte
die Eignung der Formation das CO2 über lange Zeiträume sicher verwahren zu können
Zusammenfassung iii
nachgewiesen werden. Dieser Nachweis ist nötig um den Projekterfolg sicherzustellen sowie
eventuelle Risiken zu vermeiden. Die Untersuchung dieser Fragestellungen erfolgt in dieser
Arbeit über die statistische Analyse einer Datenbank relevanter Formationsparameter,
sowie über analytische und numerische Experimente.
Ein wichtiger Aspekt ist hierbei die ausschließliche Betrachtung des Injektionszeitraumes,
d.h. es werden Zeiträume bis zu ∼20 Jahren kontinuierlicher Injektion untersucht. Dieser
Injektionszeitraum ist von hoher Bedeutung für die resultierende Speicherkapazität und
die Risikobewertung. Alle Langzeitprozesse, wie z.B. geochemische Reaktionen, hängen
von diesem anfänglichen Ausbreitungszeitraum ab, da nur in den Teilen der Formation
weitere Reaktionen stattfinden können, die auch vom CO2 erreicht werden. Ebenso ist das
Risiko für ein Austreten des CO2 aus der Formation während des Injektionszeitraumes am
höchsten, da hier der entstehende Überdruck am größten ist und sich die CO2 Fahne stetig
ausbreitet. Nach Injektionsende sinkt der entstandene Überdruck wieder ab und sekundäre
Einschlussmechanismen, wie z.B. der Einschluss des CO2 im Porenraum als nicht mehr
verdrängbare Phase aufgrund von Mehrphasen-Gesetzmäßigkeiten, die Lösung des CO2
im Formationswasser sowie die geochemische Bindung des CO2 in Mineralen, führen zu
einem beträchtlichen Sicherheitszuwachs. Zudem ist nach Abgleich gemessener Daten mit
berechneten Simulationswerten (und Anpassung der Simulationsparameter) das Vertrauen
in die Ergebnisse numerischer Prognoserechnungen hoch.
Die Entwicklung eines numerischen Modells, welches es ermöglicht die relevanten physikalis-
chen und thermodynamischen Prozesse während einer CO2 Injektion in geologische Forma-
tionen zu simulieren, war nicht Bestandteil dieser Arbeit. Es wurde ein am Institut für
Wasserbau bereits entwickeltes Modell verwendet, weshalb die Leistungsfähigkeit hier nur
kurz beschrieben wird.
Die Strömung und der Transport von CO2 in einer porösen, starren Gesteinsmatrix
werden je nach Fragestellung entweder mit einem 2-Phasen oder einem komplexeren 2-
Komponenten-2-Phasen Ansatz beschrieben. Bei beiden Ansätzen existiert eine CO2 -Phase
sowie eine Wasserphase. Im 2-Komponenten-2-Phasen Ansatz finden Lösungsprozesse
der Komponenten (also CO2 oder Wasser) in der jeweils anderen Phase statt und der
diffusive Transport wird über einen Fickschen Ansatz beschrieben. Die Beschreibung
der advektiven Flüsse erfolgt durch das Darcy-Gesetz. Die Eigenschaften der beiden
Fluidphasen sowie der eventuelle Massentransfer der Komponenten zwischen den beiden
Phasen werden in Abhängigkeit von Druck und Temperatur beschrieben. Die entste-
henden Temperaturänderungen werden im 2-Komponenten-2-Phasen Ansatz ebenfalls
berechnet. Die beschriebenen konzeptionellen Modelle, bzw. ihre mathematischen und
numerischen Realisationen, sind in das Simulationsprogramm MUFTE-UG (MUltiphase
Flow Transport and Energy Model on Unstructured Grids) implementiert. Das Modell
iv Zusammenfassung
Dimensionsanalyse
In einem ersten Schritt wurden die dominierenden Kräfte und relevanten Prozesse in einer Di-
mensionsanalyse identifiziert und bewertet. Die geltenden Mehrphasen-Bilanz-Gleichungen
wurden durch die Einführung so genannter charakteristischer Größen entdimensionalisiert.
D.h. die physikalischen Einheiten der Gleichungen werden entfernt, indem die einheiten-
behafteten Variablen relativ zu den charakteristischer Größen dargestellt werden. Um das
zu erreichen, müssen charakteristische Größen für die Länge, die Zeit, die Geschwindigkeit
und den Druck gefunden werden. Hierbei sind die Größen Geschwindigkeit, Länge und
Zeit voneinander abhängig. Die entdimensionalisierten Bilanzgleichungen wurden so umge-
formt, dass dimensionslose Kennzahlen darin zu finden sind. Diese Kennzahlen stellen
Kräfteverhältnisse im porösen Medium dar. Die Kapillarzahl stellt hierbei das Verhältnis
von Kapillarkräften zu viskosen Kräften dar und die Gravitationszahl stellt das Verhältnis
Zusammenfassung v
Ebenen eingeteilt, je nach immanenter Unsicherheit der Schätzung und nach den Kosten der
Speicherung. Auf der untersten Ebene steht die theoretische Kapazität einer Formation, dies
ist der verfügbare Porenraum abzüglich der aufgrund von Mehrphasengesetzmäßigkeiten
nicht verdrängbaren Wasseranteile. Auf der nächst höheren Ebene steht die effektive
Kapazität. Sie ist eine Untermenge der theoretischen Kapazität. Teile der theoretischen
Kapazität welche durch geologische oder ingenieurstechnische Gründe nicht erschlossen
werden können, werden hier nicht berücksichtigt. Auf den höheren Ebenen werden dann
legislatorische und Genemigungsaspekte sowie Aspekte der Infrastruktur usw. integriert.
In dieser Arbeit interessiert die Relation der theoretischen zur effektiven Kapazitätsebene.
Hier wurde ein volumenbasierter Speicherkoeffizient C definiert, der die theoretische in die
effektive Kapazität überführt. Der Speicherkoeffizient C kann noch in weitere Einzelfak-
toren zerlegt werden. Auf dieses Detail soll in dieser Zusammenfassung aber nicht weiter
eingegangen werden. In der effektiven Kapzität, welche ein Speichervolumen darstellt, kann
eine bestimmte Masse CO2 gespeichert werden, hier als Meff bezeichnet. Der Koeffizient
C und die Masse Meff wurden nun mit Hilfe numerischer 1-D und 3-D Simulationen der
definierten typischen Speicherformationen berechnet und diskutiert. Es konnten vielerlei
Aspekte der Speicherung gezeigt werden, so z.B. der Einfluss der relativen Permeabilitäten,
des Kapillardrucks oder der Injektionsrate. In der Literatur gibt es bislang keine belastbaren
Abschätzungen von C und Meff für die vorhandene Bandbreite möglicher Formationsparam-
eter, so dass diese Arbeit einen entscheidenden Fortschritt in dieser Richtung darstellt. Es
wurde weiterhin gezeigt, dass es möglich ist die berechneten Koeffizienten C und die Massen
Meff mit Hilfe der entwickelten (dimensionslosen) Kennzahlen abzuschätzen. Die in dem vo-
rangegangenen Kapitel bereits vermutete Vorteilhaftigkeit kleiner Gravitationszahlen sowie
zum Teil hoher Kapillarzahlen konnte nun quantitativ bestätigt werden. Die berechneten
Speicherkoeffizienten aller betrachteten Fälle sind kleiner als 18%, d.h. weniger als 18%
des vorhandenen Porenraums kann zur CO2 Speicherung genutzt werden. Überträgt man
den höchsten berechneten Speicherkapazitätskoeffizienten auf die Jahresproduktion eines
typischen kohlebefeuerten Kraftwerks (Ausstoß 1 MtCO2 pro Jahr) über einen Zeitraum
von 25 Jahren, so erhält man eine CO2 -Fahne mit einem Radius von 1.83 km (weitere An-
nahmen hierbei sind eine Formationsdicke von 100 m und eine CO2 -Dichte von 660.7 kg/m3 ).
Sensitivitätsanalyse
Die Frage nach dem Einfluss einzelner Formationsparameter auf das Modellergebnis wurde
in einer Sensitivtätsanalyse beantwortet. Als Modellergebnis dienten hier wieder Kriterien,
welche repräsentativ sind für die Abschätzung der Speicherkapazität einer Formation
sowie Kriterien zur Risikobewertung einer möglichen CO2 -Leckage. Die sog. Morris
Methode wurde hier angewandt um auf eine sehr effiziente Weise eine Rangfolge der
qualitativen Parametereinflüsse zu entwickeln. Diese Methode untersucht den gesamten
definierten Parameterraum, wobei jeweils ein Parameter variiert wird. Somit wird an
jedem untersuchten Punkt im Parameterraum ein lokales Sensitivitätsmaß ermittelt.
Zusammenfassung vii
Gemittelt über alle untersuchten Punkte ergibt sich daraus eine qualitative Abschätzung
des Parametereinflusses welcher in Relation zu den Einflüssen der anderen Parameter zu
bewerten ist. Betrachtet man die Standardabweichung der lokalen Sensitivitäten lässt sich,
in Relation zu den Standardabweichungen der lokalen Sensitivitäten der anderen Parameter,
eine Parameterinteraktion oder ein nicht-lineares Parameterverhalten ermitteln. Es wurden
15 Formationsparameter ausgewählt und zufällig über den gesamten, jeweils physikalisch
sinnvollen Bereich untersucht. Insgesamt wurden 64 individuelle Parameterkombinationen
erstellt und simuliert. Diese 64 Fälle zufällig generierter Parameterkombinationen bilden
die Grundlage für die Sensitivitätsanalyse. Als einflussreichste Parameter wurden die abso-
lute Permeabilität, das Injektionsintervall, und die Formationstiefe identifiziert. Hingegen
waren die CO2 Injektionstemperatur und der Fallwinkel der Formation eher vernachlässigbar.
Risikoanalyse
In einer abschließenden Risikoanalyse wurde das Risiko bezüglich einer CO2 -Leckage an
einem Brunnen in einiger Entfernung zum Injektionsbrunnen untersucht. Dabei wurden vier
einflussreiche Formationsparameter bezüglich ihres Einflusses auf das potentielle Risiko einer
Leckage statistisch untersucht.
Diese Untersuchung hatte verschiedene Ziele. Zum ersten soll sie die Möglichkeit bieten,
verschiedene Formationen in Abhängigkeit ihrer Eigenschaften bezüglich des potentiellen
Risikos miteinander zu vergleichen. Dies soll bei der Entscheidungsfindung helfen, für welche
Formationen weitere Untersuchungen vielversprechend sind. Zum zweiten soll sie bei der Po-
sitionierung des Injektionsbrunnens helfen, der von mehreren potentiellen Leckage-Brunnen
umgeben sein könnte.
Die unabhängigen Formationsparameter sind hierbei die Porosität, der geothermale Gradi-
ent, die Tiefe der Formation und die Anisotropie zwischen horizontaler und vertikaler abso-
luter Permeabilität. Diese Parameter wurden als unabhängig ausgewählt da sie einerseits in
der Sensitivitätsanlyse als einflussreiche Parameter identifiziert wurden und andererseits eine
Parameterverteilung für die statistische Untersuchung vorhanden war. Die Beschränkung
auf vier unabhängige Parameter war aufgrund des hohen Rechenaufwandes nötig den eine
Betrachtung aller einflussreichen Parameter verursacht hätte. Die weiteren Formationspa-
rameter wurden durch Funktionalitäten der unabhängigen Parameter ausgedrückt. Diese
Funktionalitäten beruhen auf einer umfassenden Literaturrecherche.
Das Risiko wurde durch die Wahrscheinlichkeit mit der eine solche Leckage auftreten könnte
berechnet, multipliziert mit der entwichenen CO2 -Masse in Abhängigkeit von der Zeit seit
Injektionsbeginn und in Abhängigkeit von der Distanz zwischen Injektionsbrunnen und
Leckage-Brunnen. Das Risiko hat somit die Einheit einer Masse. Zahlreiche Simulationen
wurden durchgeführt mit jeweils unterschiedlichen Formationseigenschaften. Die Formation-
seigenschaften wurden dabei so festgelegt, dass jeweils drei der vier unabhängigen Parameter
zufällig aus der Parameterdatenbank generiert wurden. Der vierte unabhängige Parameter,
die Anisotropie der absoluten Permeabilität, wurde zufällig aus einer Parameterverteilung
viii Zusammenfassung
generiert die auf einem theoretischen Modell basiert. Die statistischen Verteilungen der
Parameter wurden somit berücksichtigt. Die umfangreichen numerischen Simulationen wur-
den mit Hilfe des “CO2 Community Grid” durchgeführt, einer virtuellen Forschungsumge-
bung die es ermöglicht örtlich getrennte Supercomputer (Parallelrechner) in den Ländern
Dänemark, Finnland und Norwegen über ein zentrales Zugangsportal einfach und einheitlich
zu nutzen.
Aus dem berechneten Risiko wurde dann eine analytische Gleichung entwickelt. Mit dieser
Gleichung ist es möglich schnell und einfach eine quantitative Abschätzung für das durch-
schnittliche, potentielle Risiko einer Leckage zu erhalten in Abhängigkeit der Zeit und
der Distanz zum Leckage-Brunnen. Dies ist generell auch für mehrere Brunnen möglich,
sofern die Brunnen sich nicht gegenseitig beeinflussen. Des Weiteren konnten die vier un-
abhängigen Formationsparameter bezüglich ihres Einflusses auf das Risiko untersucht wer-
den. Den höchsten Einfluss hatte dabei der geothermale Gradient und die Tiefe der Forma-
tion. Ein steigender geothermaler Gradient und eine geringere Formationstiefe verursachten
ein größer werdendes Risiko zu einem gewählten Zeitpunkt, d.h. das Risiko stieg früher
an. Die Anisotropie der absoluten Permeabilität hatte einen gewissen Einfluss in einiger
Distanz des Leckage-Brunnens zum Injektionsbrunnen. Hier führte eine größer werdende
Anisotropie zu jedem gewählten Zeitpunkt zu einem größeren Risiko. Interessanterweise
hatte die Porosität keinen Einfluss auf das potentielle Risiko. Dieses Verhalten ist begründet
in der Abhängigkeit der absoluten Permeabilität von der Porosität. Mit größer werdender
Porosität, und somit erwarteter langsamerer Ausbreitung der CO2 -Fahne, verursacht durch
den vergrößerten verfügbaren Porenraum, steigt auch die absolute Permeabilität, was eine
schnellerer Ausbreitung der CO2 -Fahne vor allem in den Bereichen direkt unterhalb des
Caprocks verursacht. Diese beiden Effekte heben sich gegenseitig auf und das Risiko war
hier nahezu unabhängig von Variationen der Porosität.
Abschließend wurde eine umfassende Diskussion über die Sensitivität der getroffenen
Annahmen in Bezug auf das potentielle Risiko geführt. Dabei waren wichtige Diskus-
sionspunkte die Bedeutung der Berücksichtigung weiterer unabhängiger Parameter, die
Wahl unterschiedlicher Funktionalitäten für die Berechnung der abhängigen Parameter, die
Bedeutung der Berücksichtigung zusätzlicher Prozesse (wie z.B. geochemische Prozesse)
sowie der Einfluss realer geologischer Formationsgeometrien und Strukturen.
The severe effects that nature and mankind are facing due to increasing temperature and
changing climate include a rise in sea level, failing crop yields in many developing countries,
and the extinction of animal and plant species. This ecological and economic interest led to
the United Nations Framework Convention on Climate Change (1997), which has been ac-
cepted by 189 nations, and whose main objective is to achieve “. . . stabilisation of greenhouse
gas concentrations in the atmosphere at a level that would prevent dangerous anthropogenic
interference with the climate system. Such a level should be achieved within a time-frame
sufficient to allow ecosystems to adapt naturally to climate change . . . ”. In 2005, a num-
ber of nations approved an addition to this treaty which also includes some legally binding
measures, known as the Kyoto Protocol. Nevertheless, due to the complexity of the prob-
lem, it is not clear to date what a sustainable level of greenhouse gas concentrations is and
what emission reductions are necessary. From model runs, it has become clear that emis-
sion reductions in the range of 55–90 % by 2100, compared to the emissions of 2001, might
be necessary to stabilise the atmospheric CO2 concentration at a value of 450 ppm (IPCC,
2005). Beside the ecological interest, there is an economic interest. Stern (2007) states that
∗
parts per million, i.e. ratio of the number of molecules of the considered gas compared to the total
number of molecules of dry air.
1
2
8000
4000
2000
0
1850 1900 1950 2000
Year AD
Figure 1.1: Variation of global annual carbon emissions from burning solids (e.g. coal), liquids
(e.g. petroleum), natural gas, and others (i.e. cement production and natural gas lost
during oil and gas mining) from 1850 to 2005 (after Marland et al. (2008)). One tonne
of carbon compares to 3.6 tonnes of CO2 .
“the benefits of strong, early action on climate change outweigh the costs”. “Costs” refers
here to a global reduction of the gross domestic product.
The options for reducing global CO2 emissions are manifold, although, considering the mag-
nitude of the problem, one single option is not sufficient. Pacala and Socolow (2004) estimate
the magnitude of the problem, discuss possible options for solving it, and state that viable
techniques already exist. These authors assume that the current carbon emissions will con-
tinue to grow linearly and reach a value of 14000 Mt C/year by 2054. To stabilise the CO2
concentration at 500 ppm, it is assumed sufficient to maintain current emissions (∼7000 Mt
C/year) for the next 50 years and reduce them significantly afterwards. Thus, the total mass
of future emissions that need to be avoided in the next 50 years accumulates to 175 Gt C. The
authors then discuss 15 options for activities that reduce emissions to the atmosphere. All
necessary technologies are currently deployed on an industrial scale, but need to be upscaled.
In their concept, activities start in 2004 at zero prevented emissions and reach a value of
1000 Mt C/year prevented emissions in 2054. If seven of those 15 activities could be scaled
up to such a magnitude, this would solve the problem. Options are grouped in different sec-
tors, i.e. energy efficiency and conservation, fuel shift, carbon dioxide capture and storage,
nuclear fission, forests and agricultural soils. In conclusion, no single activity can prevent
future emissions sufficiently, but there are a number of options which can be scaled up and
simultaneously need to be expanded. The focus of this study lies on CO2 capture and storage
(CCS). The special attribute of CCS is the possibility of a fast and large-scale deployment
that could outweight time delays in the development of other technologies. The long-term
1 Introduction 3
Carbon dioxide capture and storage is a process that captures CO2 from the burning of fossil
and/or renewable fuels and from processing industries and then stores the CO2 away from
the atmosphere for geological periods of time. Carbon dioxide is primarily captured at point
sources such as power plants and other large-scale industrial processes.
To capture the CO2 , there are different approaches, i.e. pre-combustion, post-combustion,
and oxyfuel combustion. The process conditions for operation determine the approach to be
selected. Each approach requires at some point a separation of CO2 , water, or oxygen from
a bulk gas stream. Efficiencies of 80–90 % of captured CO2 can be reached whereas about
10–40 % more energy is required for the additional operations (IPCC, 2005).
It is preferable to transport the captured CO2 to feasible storage options in pipelines, but
transportation by ship, rail, or road tankers is also possible. Challenges at the transport
stage include costs, design, and safety, although experience with the current practice sug-
gests that these problems can be solved.
Carbon dioxide storage can occur in geological formations, in the ocean, or by industrial fixa-
tion in inorganic carbonates. In the following, the focus is on storage in geological formations.
Feasible reservoir types include deep saline aquifers, oil and gas fields and unminable coal
seams. Accordingly, one can distinguish between the exclusive purpose of storing the CO2
and the use of injected CO2 to enhance oil recovery (EOR), natural gas (EGR) or methane
from coal beds (ECBM). The latter three options are likely to be implemented in the near
future because of cost benefits and presumably good knowledge about the site-specific geol-
ogy as well as the existing infrastructure. However, in the following, the focus of this study
is on storage in deep saline aquifers since estimates of the available storage capacity neces-
sary to store the immense amounts of human CO2 production are promising. In Figure 1.2,
the technology is schematically sketched, including processes occurring in the sub-surface,
monitoring devices, and potential leakage pathways.
Carbon dioxide is injected into a saline formation at a depth preferably greater than 800 m
below the surface. The CO2 plume spreads laterally in the aquifer, displacing the resident
brine, which results in pressure increase. At the same time, due to the lower CO2 den-
sity than the brine density at these pressures and temperatures, it migrates in an upward
direction. To prevent CO2 from leaving the formation a confining layer above the storage
reservoir is necessary. This confining layer is usually called caprock and should provide low
permeability, considerable thickness, and no geological weaknesses such as e.g. fractures or
faults. However, the risk remains that CO2 might leak out of the storage reservoir through
these natural or man-made pathways. This issue needs to be addressed carefully in every
storage attempt. Carbon dioxide leakage is illustrated in Figure 1.2 through a leaky well
(e.g. poorly plugged abandoned well or old oil well) and through a fault. It may then mi-
4
Figure 1.2: Principal processes, leakage risks, and monitoring techniques associated with CO2
storage in geological formations (Figure courtesy of Lawrence Berkeley National Lab-
oratory).
grate into shallower aquifers (harming potable groundwater) or migrate back to the surface
(leading to risk directly associated with exposure to leaked CO2 ). To prevent harm to the
health of humans and animals as well as to the environment, an efficient and reliable moni-
toring network is essential. Monitoring leaked CO2 on the land surface, geophysical (seismic)
monitoring, and monitoring from aeroplanes is also sketched in Figure 1.2.
Currently, several projects all over the world are injecting CO2 into saline aquifers for either
socio-economic reasons or for research purposes. The first commercial attempt is being made
by the Statoil-operated Sleipner project (Torp and Gale, 2004). Since 1996, approximately
1 Mt CO2 /year has been injected into the 50–250 m-thick Utsira formation in the North Sea
at ∼1100 m depth. The CO2 is extracted from natural gas (containing about 9% CO2 ) that
is captured from another field and then processed to the supercritical conditions of 80 bar
pressure and 40 °C temperature before being re-injected. At another site in Norway, CO2 has
been injected since May 2008 - the Snøhvit field in the Barent Sea. The CO2 content of the
natural gas extracted there is decreased from 5–8 % to less than 50 ppm before the gas can
be further processed (converted to liquefied natural gas). The ∼0.75 Mt CO2 produced per
1 Introduction 5
year are re-injected into a deeper formation. Another commercial example is the In-Salah
project in the southern Sahara (Algeria), where CO2 has been injected since 2004. Similar
to the Norwegian projects, the natural gas produced initially has a CO2 content of ∼10 %,
which has to be decreased to ∼0.3 % to meet European market standards. The annually
produced ∼1.2 Mt CO2 are re-injected into a 1800 m-deep sandstone reservoir. The next
commercial project in operation may be the Gorgon Joint Venture project (Australia). The
natural gas produced there has a CO2 content of up to 14 %, which is to be reduced. The
∼2.7–3.2 Mt CO2 produced annually are to be re-injected into a saline formation at ∼2300 m
depth (International Energy Agency Greenhouse Gas R&D Programme, 2008).
Beside that, several pilot sites are currently being investigated for experimental research
purposes.
In the Nagaoka project, 10400 t CO2 were injected in 2003 and 2004 in a ∼1100 m deep
saline aquifer at the Iwanohara base near Nagaoka (Japan). The purpose of the project
was the investigation of the behaviour of CO2 during and after injection, the long-term
stability of CO2 in the reservoir, and the potential and costs of CO2 storage (Nagaoka
project consortium, 2009).
In 2004, the Frio project injected 1600 t CO2 and 320 t CO2 in two stages in two saline
aquifers at the Frio site, north-east of Houston in the U.S.A. in ∼1600 m depth. Extensive
monitoring techniques have been tested (Hovorka et al., 2006).
In the CO2SINK project, 8450 t CO2 were injected until February 8th 2009 into a saline
aquifer in the Ketzin anticline close to Berlin (Germany) at a depth of ∼500 m–700 m. It
is planned to inject up to 60000 t CO2 and sophisticated monitoring techniques are to be
tested. The project provides an in-situ laboratory for CO2 storage to fill the gap between
the numerous conceptual engineering and scientific studies on geological storage and a fully-
fledged on-shore storage demonstration (CO2SINK project consortium, 2009).
The Department of Energy in the U.S.A. has initiated a national network of seven regional
partnerships to investigate the best approaches for capturing and storing gases that can
contribute to global climate change (NETL, 2009). The partnerships aim at injecting CO2
into 14 formations in the U.S.A. in 2009 and 2010. Injection rates at the sites vary between
3 kt and 10.8 Mt CO2 .
Having outlined the motivation for this study and the current state of the CSS technology,
the principle trapping mechanisms of CO2 injected in geological formations are explained in
the following since they are essential comprehending of this study’s objectives.
very different time scales with varying contributions to the total trapping. The trapping
mechanisms are a result of fundamental processes occurring in a reservoir after CO2 has
been injected.
Common structural traps are distorted geological layers, forming e.g. folds or anticlines.
Even closed faults can act as structural traps. Stratigraphic traps are referred to if the per-
meability changes to much lower values within the respective reservoir, thus forming a seal.
The storage security, especially in the early stages after the start of the injection, is largely
influenced by the caprock integrity. A fracturing of the caprock or fault re-activation due to
over-pressurisation should be prevented by all means.
Residual Trapping
The minimum saturation that is attainable for a fluid when displaced by another (immisci-
ble) fluid from a porous medium is called the residual saturation. The residual saturation
cannot be further reduced by viscous forces. Thus, the fraction of the CO2 plume below
residual saturation is residually trapped. This is of special importance in regions where the
CO2 plume is displaced again by brine. The CO2 in residual saturation may then dissolve in
brine and eventually CO2 as a separate phase disappears. The effect of hysteresis (discussed
in Section 2.2.4) is also of importance here, since residual saturations may change, possibly
increase, with every drainage-imbibition cycle. Such a drainage-imbibition cycle could be
induced by non-continuous or varying CO2 injection rates.
Solubility Trapping
When CO2 comes into contact with resident brine, it immediately begins to dissolve up to
a solubility limit. Thus, considerable amounts of CO2 dissolve over time and by contact
with fresh brine as the CO2 plume expands. Since a de-pressurisation of the brine with CO2
load is not to be expected, CO2 is safely stored. This mechanism is referred to as solubility
trapping. Since brine with dissolved CO2 has a higher density than fresh brine, it slowly
migrates to deeper regions and fresh brine may be available for further dissolution.
Mineral Trapping
When CO2 is dissolved in brine, it may form ionic species as the rock dissolves, accom-
panied by an increase in pH. Finally, some fraction may be converted to stable carbonate
minerals. This process is called mineral trapping (IPCC, 2005). Significant amounts of CO2
1 Introduction 7
mass trapped in minerals is expected only after some hundreds to thousands of years. How-
ever, chemical reactions are not included in this study due to the focus on the injection stage.
Time Scales
On a short time scale, structural and stratigraphic trapping are the dominant trapping mech-
anisms. Over time, the contribution of residual, solubility and mineral-trapping mechanisms
increases, as does the storage security, since these mechanisms represent a long-term fixation
of CO2 (IPCC, 2005). The mineral-trapping mechanism is presumably negligible here since
its contribution to trapping may only become significant after centuries. The different im-
portance of the different trapping mechanisms over time is due to changes of the dominant
processes. The advection-dominated multiphase flow regime occurs on the short to medium
time scale. With declining pressure gradients and density differences, these processes lose
their driving forces. Inter-phase mass-transfer processes, like the dissolution of CO2 in brine
gain importance on a large time scale. Geochemical reactions eventually lead to a very secure
mineral trapping of CO2 . Figure 1.3 schematically shows the contribution of the individ-
ual trapping mechanisms to the total trapping versus time and variation of the dominant
processes.
100
trapping contribution %
residual
trapping
increasing storage security
solubility
trapping
mineral
trapping
0
geochemical
reactions
phase transfer processes
multiphase
behavior dissolution and diffusion
advection−dominated
(viscous, buoyant, capillary)
Figure 1.3: Storage security depending on different trapping mechanisms and dominant processes
over time (modified after IPCC (2005), Class (2008)).
8 1.2 Objective of the Study
… Proposal to regulator •…
… Site closure •…
Figure 1.4: Life cycle of a CO2 storage project (modified after IPCC (2005)).
This study refers to the key questions in the early phase of site screening and characteri-
sation, i.e. “Is storage capacity likely to be adequate?” and “Is the site suitable?”. Site
screening is the initial step in site characterisation and is aimed towards a pre-selection of
potential storage reservoirs. Usually, little information is available on the reservoir proper-
ties and geology at this stage of a project. Consequently, the entire range of parameters
is considered in this study. At a later stage, further investigations on the properties of the
identified reservoirs will lead to good data availability, which can then be fed into detailed
investigation methods, e.g. site-specific numerical models. This study does not aim at devel-
oping a comprehensive framework which can be applied by regulators or operators. Rather,
fundamental research is conducted and the underlying processes of CO2 injection in geolog-
ical formations are investigated. This is done by keeping the focus on the key questions,
i.e. estimating storage capacity and evaluating risk. In site screening, operators usually try
to identify sites with an anticlinal shape, high porosity, and high permeability, at a depth
1 Introduction 9
of around 1000–2000 m. Such a site is potentially of low risk (since CO2 spreading is lim-
ited by the anticlinal structure), provides good injectivity to inject CO2 economically, and
the estimated storage capacity is high. However, these sites are limited in number and will
therefore not be available for a large-scale deployment of this technology after some time.
Hence, the aquifers investigated here have a horizontal caprock and are limited by a defined
leakage point at some distance to the injection point.
Following this argumentation, the outline of this study is to introduce the conceptual, math-
ematical, and numerical model to investigate the processes of interest and then to inves-
tigate the statistical characteristics of the relevant parameter distributions as a basis for
analytical and numerical investigations. A dimensional analysis of the governing equations
is subsequently conducted to identify and assess the dominant forces and processes. The
carbon dioxide storage capacity is then defined and evaluated, based on the identified typ-
ical reservoir properties. It is analysed by using the dimensionless indicators developed. In
a sensitivity analysis the influence of various reservoir parameters on the model results is
assessed with respect to storage capacity and risk. Subsequently, a risk analysis investigates
the effect of the parameters with the greatest influence on the risk of CO2 leakage through
a pre-existing well. A risk surface is calculated to estimate average risk quickly depending
on the distance between the leaky well and the injection well and time.
The time scale considered here covers the injection stage. Processes in the early stages are
presumably of higher relevance to storage capacity. All long-term processes, like geochemi-
cal reactions etc., depend on these early events since they only occur where CO2 is present.
Similarly, the risk of leakage is highest during the injection stage due to increasing reser-
voir over-pressure and growing CO2 plume size. After shut-in, pressure recovers, secondary
CO2 trapping mechanisms contribute significantly to storing the injected CO2 safely and
confidence in history-matched numerical models is high (see Figure 1.5† ). One exception to
this reasoning is geochemical reactions, which could reduce caprock sealing efficiency. These
processes can lead to increased risk in the long term.
Figure 1.5: Hypothetical risk profile during the life cycle of a CCS project. Risk of leakage in-
creases due to increasing reservoir over-pressure and growing CO2 plume size. After
shut-in, pressure recovers and secondary CO2 trapping mechanisms contribute signif-
icantly to storing injected CO2 safely.
storage in saline aquifers has been simulated. Numerical simulation is used to answer mani-
fold questions related to the description of the processes occurring in the reservoir, e.g. stor-
age capacity estimation, design of the injection site, design of the injection strategy, design of
the monitoring network, assisting risk assessment by quantification of possible leakage rates,
estimation of over-pressures, etc. Current fields of research include the coupling of multi-
phase flow and transport with other simulators describing geomechanical and geochemical
processes and the proper description of the fluid properties of complex gas mixtures (e.g. CO2
with some impurities like H2 S) (International Energy Agency Greenhouse Gas R&D Pro-
gramme, 2008). Code comparison studies have been conducted to test consistency between
different approaches (Pruess et al., 2003).
Geochemical processes
Numerical models have been developed to describe chemical reactions occurring in CO2
storage by e.g. Gunter et al. (1997), Xu and Pruess (1998), and Clauser (2003). Questions
1 Introduction 11
of interest include long term storage capacity contribution by mineral trapping (Kühn and
Clauser (2006), Xu et al. (2007), Mito et al. (2008)), caprock alterations as a result of chem-
ical reactions (Gherardi et al., 2007), and possible injectivity reduction. Fields of current
research include the collection of basic thermodynamic and kinetic data at elevated pressure
and temperature conditions and the modelling of reactions assuming impurities in the CO2
stream (Gaus et al., 2008).
Geomechanical processes
The description of geomechanical processes occurring as a result of the injection of CO2 into
saline aquifers is undertaken by coupling geomechanical models to flow and transport models
(Le Gallo et al. (2006), Rutqvist et al. (2008)). Various approaches have been followed for the
coupling strategy, either a close coupling (feedback of geomechanical effects on porosity and
permeability and thus on the flow and transport processes) or a loose coupling (no feedback
effects and sole description of stresses and possible effects thereof). Questions of interest in
geomechanical models are concerned with the reactivation of existing faults and hence the
estimation of a maximum sustainable injection pressure, and the estimation of porosity and
permeability changes and hence the influence on the flow and transport processes.
Dimensional Analysis
Dimensional analysis is a mathematical procedure for describing complex system behaviour
in an elegant way. All physically meaningful equations can be converted into a dimen-
12 1.3 State of the Art
Storage-capacity assessment
As outlined by Bachu et al. (2007), a variety of approaches and methodologies for assessing
CO2 storage capacity has led to conflicting results for local, regional and global estimates.
At present, two major methodologies are proposed. Bachu et al. (2007) propose the concept
of so-called “Resource-Reserve Pyramids” in which several aspects of CO2 storage are con-
sidered, e.g. time scales, assessment scales, and various storage options (in saline aquifers,
in enhanced gas recovery, etc.). For example the “Techno-Economic Resource-Reserve”
pyramid consists of theoretical, effective, practical, and matched storage-capacity estimates,
having the units of a volume. The uncertainty of the estimates is reflected by its place in the
pyramid. The theoretical capacity is the entire pore space of a formation, reduced by the
residual (irreducible) water fraction. Effective storage capacity is the estimate of interest in
this study, since it is a subset of the theoretical capacity estimate, which satisfies a range of
geological and engineering constraints and which can be estimated by numerical simulation.
This corresponds to the term “resources”. Computationally, the effective capacity can be
estimated by multiplying the theoretical capacity by a capacity coefficient. Estimations of
this capacity coefficient are, however, not given by the authors. The practical storage ca-
pacity is defined as the reserves, considering economic, legal and regulatory, infrastructure,
and general economic aspects. The matched capacity then results from detailed source–sink
matching. To obtain an estimate related to the mass of CO2 that can be stored, the volume
estimates are multiplied by the average in-situ CO2 density. The concept proposed by NETL
(2007) is computationally equivalent to the effective capacity estimate of Bachu et al. (2007)
for saline formations if the residual (irreducible) water fraction is included in the capacity
coefficient and if the average CO2 density is used in the equation rather than minimum and
maximum values (International Energy Agency Greenhouse Gas R&D Programme, 2008).
As an advancement, the NETL (2007) approach provides estimates for the capacity coeffi-
cient. The estimates were gained by assuming minimum and maximum values for the relevant
reservoir properties and processes (e.g. porosity, height of the reservoir, areal displacement
efficiency, vertical displacement efficiency, etc.), with normal distributions of these proper-
1 Introduction 13
ties and after Monte Carlo simulations (no numerical multi-phase flow simulations). The
resulting capacity coefficient ranges between 1 and 4 % for the 15 – 85 % confidence interval.
As an assumption, the CO2 injection wells are placed regularly throughout the basin/region
to maximise storage.
Sensitivity analysis
Sensitivity analysis in general is the study of how a variation in the model input translates
into a variation in the model output. As such, it is a method that is widely applied in various
fields. Saltelli et al. (2000) review various methodologies and evaluate the trade-off between
the computational cost and the possible information to be gained. In this study, the method
proposed by Morris (1991) and extensions thereof are applied. Campolongo et al. (2004)
tested the Morris Method and stated that the method is efficient in identifying irrelevant
model input factors.
Risk analysis
Risk analysis is in general the systematic study of how uncertainties propagate into risks
encountered and the estimation of the impact. It assists in developing mitigation options
for minimising or preventing harm. With respect to CO2 storage, there is considerable un-
certainty involved in the knowledge about the subsurface. Consequently, the risk of CO2
leaking out of the storage formation is a major concern which needs to be taken care of.
Celia and Bachu (2003) state that in North America, due to the large number of existing
wells and due to the sparse knowledge about the state of those wells, the leakage of CO2
does not seem to be avoidable, given the large spatial and temporal scales associated with
CCS. Subsequently, Gasda et al. (2004) analysed the spatial characteristics of existing well
locations in the Alberta basin (Canada) and state that a typical CO2 plume can encounter
up to several hundreds of wells in high-density areas. Nordbotten et al. (2004) developed
the already mentioned semi-analytical solutions to efficiently describe plume evolution and
leakage in these complex aquifer systems. This yields a useful tool for estimating risk for
systems where risk is mainly associated with the existing wells and their largely unknown
state. A more general and comprehensive screening and ranking framework is presented by
Oldenburg (2007) for evaluating potential storage sites on the basis of risk to health, safety
and the environment arising from possible CO2 leakage. The framework is based on expert
judgement on e.g. the relative importance, relative risk, and certainty of an assessment at-
tribute. However, up to now, there is no consistent risk-assessment methodology for CCS
projects.
and the outcome just serves as a guideline. International Energy Agency Greenhouse Gas
R&D Programme (2008) extended this set by indicators like injectivity, site logistics (e.g. dis-
tance from CO2 source), and potentially compromised natural resources (e.g. groundwater,
coal, etc.). Site screening can be seen as the initial step to site characterisation, leading
ultimately to site selection. After Friedmann (2007) it includes the characterisation of in-
jectivity, storage capacity, and effectiveness (long-term ability to store CO2 ). It involves
detailed geological characterisation, numerical flow and transport modelling, geochemical
and geomechanical assessment, risk assessment, monitoring programme design, and trans-
port (International Energy Agency Greenhouse Gas R&D Programme, 2008). The approach
by Oldenburg (2007) already mentioned may be considered as a site screening tool, based on
risk. Site screening and in particular site characterisation depend heavily on data availability
(e.g. maps, seismic data, well logs, etc.).
2 Conceptual, Mathematical and
Numerical Model
One important step in numerical simulation is the development of a conceptual model, the
translation of the conceptual model into a mathematical model, and the implementation
into a numerical model. The conceptual model describes the essential features and principal
processes of the system of interest, but simplifies or neglects minor features and processes
so that the resulting computational cost of the simulation can be handled. The degree of
abstraction of the real world system is therefore dependent on the scope of the investigation,
but also on other constraints, like the level of knowledge about the system, the spatial and
temporal scale of interest, and the availability of computational power. Nevertheless, it is
most important not to over-simplify the system and to guarantee the essential features and
principal processes. As previously stated, the focus of this study is the investigation of
CO2 injection processes in geological formations for site screening. The temporal scale of
interest covers the injection stage. The advection dominated multi-phase flow and transport
processes are of importance on this scale, together with diffusion, and dissolution of CO2 into
the resident brine. In site screening, it is assumed that only very few data about the geological
formations are available. Hence, chemical reactions, contributing significant amounts to the
overall CO2 trapping only after some decades and requiring knowledge about the mineral
composition of the reservoir rock, are not considered. Furthermore, geo-mechanical processes
are neglected, since it is assumed that appropriate measures to prevent caprock fracturing
can be applied when designing the injection facility. For such a system the conceptual model
is defined in the following. The mathematical and numerical model is briefly introduced.
The interested reader is however referred to Appendix A, since a detailed discussion may
detract the readers attention from the focus of this study (and the model development was
not part of this study).
15
16 2.1 Basic Definitions
this interface, mass transfer can occur from one phase to the other. Generally, one can
distinguish between solid, liquid and gaseous phases. Several solid or liquid phases may
coexist, whereas only one gas phase may occur. Each phase can be composed of one or
several components. A component relates here to either a pure chemical substance or a
mixture of chemical substances. To consider a mixture of substances, different components
are lumped into a so-called pseudo-component. A pseudo-component employs (fraction-
weighted) average properties of the pure substances and considers the interaction effects
between the substances as well.
Any chemical substance can occur in a solid, liquid, and gaseous state of aggregation. A
transition from one state of aggregation to another is possible by change of the so-called
state variables. State variables are, for instance, temperature, pressure, mass, volume. The
state variables characterise the thermodynamic state of the system and are independent
of any previous state of the system. Generally, one can distinguish between intensive and
extensive state variables. Intensive state variables are independent of the size of the system,
hence temperature and pressure are intensive state variables. Extensive state variables,
however, vary with system size, e.g. mass and volume. When dividing two extensive state
variables, a specific intensive state variable is the result, e.g. specific volume= volume
mass
. If
the state of aggregation changes from one state to another, a phase-coexistence curve (also
known as phase boundary) is crossed. At this curve, defined by a specific set of state
variable values, the substance forms two stable phases. One example is the phase-coexistence
curve that separates the gaseous from the liquid state of aggregation, the so-called vapour-
pressure curve. When crossing the phase-coexistence curves, fluid properties usually change
discontinuously. In Figure 2.1, a schematic diagram shows the different states of aggregation
for CO2 .
Three phases coexist at the end points of a phase-coexistence curve. In Figure 2.1, this
occurs at the critical point and at the triple point. At temperatures and pressures above
the critical values (pcrit,CO2 = 7.38 MPa = 73.8 bar, Tcrit,CO2 = 304.1 K = 30.95 ◦ C), the fluid
is in a supercritical state of aggregation. This state of aggregation is of special interest for
CO2 in storage scenarios, since the resulting high fluid density allows an efficient utilisation
of the available pore volume.
In a multi-phase system, the states of aggregation of the individual phases may change in-
dependently. Hence, different combinations of states of aggregation for the phases may occur.
2 Conceptual, Mathematical and Numerical Model 17
p
supercritical
[bar]
73.8
(pcrit) critical
liquid point
solid
5.11
(p3 ) triple− gaseous
point
−56.35 30.95
(T3 ) (Tcrit ) T [°C]
Figure 2.1: Schematic phase diagram of carbon dioxide, depending on pressure (p) and temper-
ature (T). The liquid-gas phase boundary between the triple point and the critical
point is known as the vapour-pressure curve (Lüdecke and Lüdecke, 2000).
In a multi-phase system, it is possible that fluid phases appear or disappear, e.g. due to
displacement or due to mass transfer processes. A change in the number of locally present
phases present is called a phase change. In the system investigated here, two fluid phases
are considered, CO2 and water. In each point in space, both fluid phases may occur, or only
one of them.
In the second, more complex model (cf. Section A.2), one solid phase, one liquid water-
component-rich phase (called water-rich phase in the following) and one either liquid or
gaseous CO2 -component-rich phase (called CO2 -rich phase in the following) are considered.
The water-rich phase can consist of the pseudo-component brine and a pure CO2 component.
The CO2 -rich phase consists of a CO2 component only.
18 2.1 Basic Definitions
The equation of state (EOS) for a pure substance is a mathematical formulation describing
the equilibrium relationship between pressure, temperature, and volume. Figure 2.2 shows a
schematic diagram of such a relationship. Figure 2.1 is a projection of the pVT-surface onto
the pT-plane. Van Der Waals (1873) proposed the first EOS as a polynomial of third order.
Peng and Robinson (1976) extended this approach by modifying the original equation and
introducing additional parameters. By parameter fitting, a good match can be obtained for
most substances between the EOS given by Peng and Robinson (1976) and experimental
results.
solid−liquid
sol. sup
erc
riti
CP cal
liq.
supercritical gas
p liq.+ga eou
CP T s. s
1
0
supercritical
liq. solid+g
aseous
solid gaseous Tme
Tb CP V
solid
p liq.
1
0
liq.+ga
s.
p gaseous
T 1 0
0
1 1 2 Tb CP : critical point
1
0 Ts
3 1−2−3 : triple line
solid+ 2−CP : boiling line
gaseou
s
V Ttrip 3−CP : dew line
Tme : melting temperature
Tb : boiling temperature
T3 : triple point
liq. CP
supercritical
solid−liquid liq.+ga
s.
gaseous
solid
+gase
ous
solid V T
Figure 2.2: Schematic diagram of the equilibrium relationship between pressure (p), volume (V),
and temperature (T) (known as the pVT-surface) of a pure substance which contracts
upon freezing and projections onto the pT-, pV-, and VT-diagrams (Bielinski, 2006).
Mixtures of substances can also be described by an EOS by setting up mixing rules, taking
into account the properties of the pure substances and the interaction effects between them.
2 Conceptual, Mathematical and Numerical Model 19
nC
α
xC
α = P C
, (2.1)
C nα
where nC
α is the number of moles of component C in phase α. Similarly, the mass fraction
X of component C in phase α is defined as
mC
XαC = P α C , (2.2)
C mα
where mC
α is the mass of component C in phase α.
Mole and mass fractions add up to unity for each phase by definition:
X X
xC
α = XαC = 1. (2.3)
C C
The molecular weight M c is the relation between these mass and mole fractions:
M c = mC C
α /nα . (2.4)
In the following, only mass fractions are used. In this study, the mass fractions of CO2 in
the water-rich phase XwCO2 , and the mass fraction of water in the CO2 -rich phase XCO2
w
are
of importance.
2.1.6 Salinity
In deep geological formations, considerable amounts of salt are dissolved in the formation
water. It is therefore often referred to as brine. To describe the amount of salt dissolved in
the water-rich phase, the term “salinity” is introduced. Different definitions of salinity exist.
In this study, salinity is defined as the mass fraction of salt related to the total mass of the
solution [kg salt/kg solution]:
msalt
S= , (2.5)
msalt + mw
km m mm µm
Figure 2.3: Different scales for fluid flow in porous media (modified after Kobus et al. (1996)).
On the molecular scale [∼ 10−10 m] the movement of individual molecules and the inter-
action with other molecules can be described. When such a system is modelled (e.g. in
biochemistry), the computational cost is extremely high, even for small systems due
to the large number of molecules (6.0221 · 1023 ) per mole of any substance.
On the microscale [∼ 10−3 m] an averaging is applied over individual molecules and their
interactions. The “continuum approach” describes this averaging procedure over a
sufficiently large number of molecules, consequently they can be assumed to be con-
tinuously distributed in space. Thus, new variables appear, e.g. density and viscosity.
By solving the Navier-Stokes equations, multi-phase fluid flow can be described over
a volume of several pores. Due to a lack of knowledge about realistic pore structures
or the high cost of implementing these structures, simplified pore-network models are
often used.
On the macroscale [∼ 10−1 m] the microscopic properties of the system are averaged within
a defined volume. This volume is called the “representative elementary volume (REV)”.
It has to meet the requirement that the averaged properties are independent of minor
2 Conceptual, Mathematical and Numerical Model 21
changes in volume size. Figure 2.4 shows schematically the definition of a REV for the
volume fraction of the pore space within the respective volume, known as porosity on
the macroscale. For small volume sizes, the pore space fraction fluctuates. When the
volume is increased up to a certain minimum REV size, indicated by V0 , the pore space
fraction (porosity) remains constant. When the volume is further increased, there is a
limit where heterogeneities start to influence the averaging procedure; this represents
the maximum REV size. As already indicated, new variables appear as a result of the
averaging on the macroscale, these are e.g. porosity and saturations of the individual
phases. This means that the distribution of the fluids in the pore space is not described
exactly any more, but expressed by the volume fraction Sα the fluid occupies within
the pore space of the respective volume. This is expressed in Equation 2.6:
heterogeneous
medium
Porosity
homogeneous
medium
0
V0=REV Volume V
Figure 2.4: Definition of a representative elementary volume modified after Bear (1972).
On the field scale [∼ 102 m] no additional averaging procedure occurs. The properties of
this scale result from the averaging procedure described previously on the macroscale
and the occurring heterogeneities are described by the REV properties. The field scale
describes the scale of interest in this study. A typical geological formation considered
for CO2 storage extends laterally up to several kilometers. In a vertical direction, the
scale of interest ranges from a few meters (for a thin injection formation only) up to
some hundred meters (by considering several formations).
22 2.2 System Properties
p V = n R T, (2.10)
where R is the ideal gas constant (8.31447 J/(mol K)). With the identities given in Equa-
tions 2.4 and 2.9, the ideal gas law can be reformulated to
pM
%mass = . (2.11)
RT
Thus, the density of an ideal gas increases with increasing pressure and decreasing tempera-
ture. The ideal gas law represents an equation of state (cf. Section 2.1.4) which neglects the
size of the molecules and the intermolecular attractions. Thus, it is most accurate at at low
pressures (i.e. large volumes) and at high temperatures (i.e. high thermal kinetic energy).
The density of carbon dioxide and water at the conditions of interest here, i.e. high pressures
and temperatures, cannot be described by the ideal gas law. In the following, the approaches
for calculating the phase densities are given as used in the model.
Density [kg/m3]
0 200 400 600 800 1000
0 0
Density
Dynamic Viscosity
500 50
Depth below Surface [m]
1000 100
Pressure [bar] *
150
1500
* Pressure is
hydrostatic
with a const. 200
2000 brine density 3
of 1100 kg/m
Geothermal 250
2500 Gradients:
High 0.062 °C/m
Median 0.030 °C/m 300
Low 0.018 °C/m
3000
High Median Low
G. Grad. G. Grad. G. Grad. 350
3500
0 20 40 60 80 100
Dynamic Viscosity [MPa s]
Figure 2.5: Variation of carbon-dioxide density and dynamic viscosity with depth for various
geothermal gradients (high geothermal gradient = 0.062 ◦ C/m, median geothermal
gradient = 0.03 ◦ C/m, low geothermal gradient = 0.018 ◦ C/m). A surface temper-
ature of 10 ◦ C and a hydrostatic pressure distribution corresponding to a water-rich
phase with density of 1100 kg/m3 are assumed.
Figure 2.6 shows variations in the water density as a result of varying pressure, temperature,
and salinity. Density decreases with pressure and salinity, and increases with temperature.
24 2.2 System Properties
Figure 2.6: Variation of water density with pressure p and temperature T for salinity values of
0.0 (bottom plane), 0.15 (centre plane), and 0.3 (top plane).
Viscosity
Viscosity is a fluid’s internal resistance to flow. The dynamic viscosity µ is the relation of the
fluid’s shear stress τ and the velocity gradient in the direction perpendicular to the direction
of flow dv/dn
τ
µ= . (2.13)
dv/dn
For some cases, it may be advantageous to use the kinematic viscosity ν, which is dynamic
viscosity divided by the fluid’s density (ν = µ/%). Generally, the viscosity of a fluid is a
function of pressure, temperature, and its composition.
-2.5
Log 10 (Viscosity/1 Pa s)
Log10 (Viscosity) [Pa s]
-3
-3.5
0
0
25
0.1 50
C]
S [-] 0.2 75 T [°
0.3 100
Figure 2.7: Variation of water dynamic viscosity with temperature T and salinity S.
Enthalpy
To describe non-isothermal processes occurring during the injection of CO2 into porous me-
dia, the caloric state variables enthalpy and internal energy are introduced. Internal energy
is related to the molecular structure and the degree of molecular activity in a system. For a
gaseous fluid, the internal energy is mainly composed of kinetic energy due to motion of the
molecules. For a solid or liquid fluid, considerable contributions are also made by the poten-
tial energy of attraction or repulsion in between individual molecules. The extensive state
variable enthalpy (H) can be described by adding volume-changing work to the extensive
internal energy (U ):
H = U + p V. (2.14)
and density variations with pressure are high. For liquids having only a low compressibility,
the influence may be of minor importance. For mixtures of substances, the individual con-
tributions of the components need to be considered as well as the heat of dissolution of the
components.
100
50
Specific Enthalpy [kJ/kg]
-50
-100 80°C
-150 60°C
-200 40°C
30°C
-250 20°C
10°C
-300
50 100 150 200
Pressure [bar]
Figure 2.8: Variation of specific enthalpy of CO2 with pressure for various temperature levels
(after Bielinski (2006)).
The reference states of the CO2 -rich phase and the water-rich phase have to be identical
to quantify enthalpy changes correctly. Span and Wagner (1996) use a reference state of
pCO2 = 1.013 bar and T = 25 ◦ C, resulting in a constant difference in specific enthalpy of
21.91 kJ/kg. However, differences in specific enthalpy are not affected by this selection of the
reference state (Bielinski, 2006). If the vapour-pressure curve is crossed, e.g. for the 10 ◦ C
curve at a pressure of 45.5 bar, there is a discontinuous decrease in specific enthalpy. For
supercritical temperatures (or pressures), this is not observed.
2 Conceptual, Mathematical and Numerical Model 27
where hpure water is the enthalpy of pure water (dependent on pressure and temperature),
hNaCl is the enthalpy of salt (dependent on temperature), and hL indicates heat of dissolu-
tion (dependent on temperature) (Garcia, 2003). Figure 2.9 shows the variation of specific
enthalpy of water, salt, and CO2 and variation of the heat of dissolution of CO2 in the
water-rich phase with temperature at a constant pressure of 100 bar. The specific enthalpy
of pure water is calculated after IAPWS (1997) and the specific enthalpy of CO2 follows the
approach of Span and Wagner (1996). The specific enthalpy of NaCl is calculated by
Z
h = u = c(T ) dT, (2.17)
where c is the specific heat capacity of salt, taken from Daubert and Danner (1989). The
volume-changing work is neglected here, i.e salt is assumed to be incompressible.
The heat of dissolution of CO2 in water ∆hL (T )CO2 is calculated after Duan and Sun (2003)
and is negative in the temperature range of interest; thus, the dissolution causes a warming of
the water-rich phase (exothermic reaction) and overall water-rich phase enthalpy is reduced.
The heat of dissolution of salt in water ∆hL (T )NaCl (after Michaelides (1981), including the
corrections after Gudmundsson and Thrainsson (1989)) is positive for the entire tempera-
ture range considered here; hence, the dissolution causes a cooling of the water (endothermic
reaction).
To summarise, Table 2.1 gives on overview of the fluid-property dependencies and literature
sources used.
28 2.2 System Properties
200
p=100bar
hNaCl
0
-200
hCO2
-400 ∆hL,CO2
0 50 100
T [°C]
Figure 2.9: Variation of specific enthalpy and heat of dissolution of CO2 in the water-rich phase
with temperature at a constant pressure of 100 bar.
Table 2.1: Fluid properties of CO2 -rich- and water-rich phase and dependencies on temperature
(T ), pressure (p), salinity (S), and mass fraction of CO2 in the water-rich phase (XwCO2 )
(Bielinski, 2006).
2 Conceptual, Mathematical and Numerical Model 29
Porosity
Porosity has already been introduced exemplarily as a soil property on the macroscale. Total
porosity describes the volume fraction of the pore space compared to the entire volume of
the representative elementary volume (REV):
Vpore space
φ= . (2.18)
V
This study always refers to effective porosity, i.e. the fraction of total porosity that is accessi-
ble to fluid flow. In comparison to the total porosity (as described above), effective porosity
excludes e.g. pore fractions occupied by water bound to soil particles or isolated pores that
are not connected.
Permeability
The absolute (intrinsic) permeability k is a measure of the resistance of a porous medium to
transmitting fluids. Generally, the absolute permeability is dependent on the properties of
the pore space, such as the porosity and several structural parameters of the porous medium
(e.g. tortuosity of the flow channels). However, there is no exact method of calculating
permeability and in practice it needs to be measured, e.g. by pumping tests (using Darcy’s
law, cf. Section 2.3.1), or to be estimated using empirically derived dependencies (e.g. on
grain-size distribution). Approaches to estimating the absolute permeability, e.g. depending
on the pore space geometry, have been made by Pape et al. (1999). Absolute permeability
may vary over several orders of magnitude, even for a single rock type such as sandstone
(Clauser, 1992). The related hydraulic permeability kf is dependent on both the properties
of the porous medium and the properties of the fluid of interest:
%g
kf = k , (2.19)
µ
where g is gravity. Absolute and hydraulic permeability are tensorial quantities and in the
typical geological formations of interest here, they may also vary by orders of magnitude
depending on the direction of flow.
Heat Capacity
Heat capacity is a measure of the heat energy required to increase the temperature of a
substance or system. It is an extensive state variable. With respect to this study, specific
30 2.2 System Properties
heat capacity is more relevant. It is defined as the heat energy required to increase the
temperature of a unit mass of a substance or system by a certain temperature interval. Here it
is always referred to kJ/(kg K). For substances or systems with considerable compressibility,
it may be necessary to distinguish between specific heat capacity at constant pressure or at
constant volume. In this study, however, the porous medium is assumed to be incompressible,
which leads to identical specific heat capacity at constant pressure and at constant volume,
referred to as the specific heat capacity of the soil grains cs . Moreover, the energy content
of the fluid phases is described using the specific internal energy (u) (see Section 2.2.2).
Residual Saturation
The residual saturation defines the minimum saturation that is attainable for a fluid when
displaced from a porous medium by another immiscible fluid. It is a measure of the fluid
which cannot be displaced, due to viscous forces. On the microscale, the wetting phase (in
this study, always the water-rich phase) residual saturation is caused by strong capillary
forces that prevent further fluid displacement. For the non-wetting phase (in this study,
always the CO2 -rich phase), this is caused by entrapped fluid bubbles in larger pores. It is
also possible that a fluid phase is trapped in several pores if it is entirely surrounded by the
other phase. The residual saturation is not only dependent on the properties of the porous
medium and the fluids, but also on the history of displacement processes that occurred in the
system. This is discussed further in Section 2.2.4 (Hysteresis). If the fluids are considered
to be miscible, however, saturations smaller than the residual saturation can be attained for
the displaced fluid (e.g. due to dissolution or evaporation).
Capillary Pressure
Figure 2.10: Interfacial tension between CO2 and water dependent on pressure and temperature
after an equation given by Kvamme et al. (2007). Accuracy is reported to be greater
than 95 % within the experimental range. The experimental range includes temper-
atures between 278 K and 335 K and pressures between 0.1 MPa and 20 MPa. The
range included in this study (cf. Chapter 7) is expanded beyond the experimen-
tal range; however, interfacial tension does not vary significantly at high pressures
and temperatures (corresponding to great depth). At shallower depth, pressure and
temperature conditions are covered by the experimental range.
In Figure 2.11, a sketch of the effect of interfacial tensions in a porous medium is shown.
The contact angle Θ is dependent on the interfacial tensions after Youngs’s equation (Helmig,
1997). Contact angles smaller than 90 degrees indicate the wetting fluid. At the equilibrium
level h, the surface tension σ of the free liquid surface causes capillary forces Fcap that are
in equilibrium with the gravitational forces Fgrav of the wetting-phase (w) column,
2π r σ cos(Θ) = π r2 h g %w , (2.20)
where r is the radius of the capillary tube. Hence, the equilibrium level h can be formulated
as
2 σ cos(Θ)
h= . (2.21)
r g %w
p
Reformulating Equation 2.21 (h= %g
+ z) yields the capillary pressure pc
2 σ cos(Θ)
pc = . (2.22)
r
32 2.2 System Properties
2r
R
Fcap
σsw Θ Fgrav
σwn h
wetting fluid
Figure 2.11: Single capillary tube with interfacial tensions between solid and non-wetting fluid
σsn , between solid and wetting fluid σsw , and between wetting and non-wetting fluid
σwn . Capillary forces Fcap cause a rise of the wetting fluid in the capillary tube
until equilibrium is reached with gravitational forces Fgrav , indicated here at level
h. Between the fluid interface and the solid wall, the so-called contact angle Θ is
defined (after Ochs (2006)).
2
10 1
M Brooks & Corey Model M Brooks & Corey Model
V Viking Formation V Viking Formation
E Ellerslie Formation E Ellerslie Formation
101 0.8 B Basal Formation
B Basal Formation M V E
0.6
B
10-1 B 0.4 M
E E
V V B
10-2 M 0.2
10-3 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
S W [-] S CO2 [-]
Figure 2.12: Capillary pressure-saturation relations (left) after Bennion and Bachu (2006) mea-
sured for CO2 -brine systems in cores of the Alberta basin in Canada (Viking,
Ellerslie, and Basal Formations) and for a Brooks and Corey model (Brooks and
Corey, 1964) employing input parameter settings as given in Table 2.2. Relative
permeability-saturation relations (right) modified after Bennion and Bachu (2008)
for the same reservoirs described above and for the same Brooks and Corey model.
Letter coding according to Section 3.4, where these relations are used to define
typical reservoir cases (full lines: brine; broken lines: CO2 ).
Relative Permeability
Relative permeability accounts for the fact, that, if multiple fluids share the same pore space
of a porous medium to flow, they interfere with each other. This can be explained by reduced
pore throats available to each fluid to flow and by the increased tortuosity of the flow paths
(since e.g. small pores might be fully saturated by the wetting-phase and the non-wetting
phase has to circumvent these pores). However, due to the complexity of the pore geome-
try, these microscale considerations can only be represented qualitatively on the macroscale
(Helmig, 1997).
On the macroscale, the effect of relative permeability is accounted for by multiplying the
absolute permeability k by a dimensionless value kr,α , called the relative permeability, that
is specific for each phase. Since kr,α ranges between zero and one, k is always reduced. As
34 2.2 System Properties
Table 2.2: Residual saturations for water (Sw,r ) and CO2 (SCO2 ,r ) and Brooks and Corey (1964)
model parameters (λBC and pd ) for measured relative permeability-saturation relations
(Bennion and Bachu, 2008) and for a synthetic model. The Brooks and Corey model
parameters (λBC ) for the measured relative permeability relations were obtained by
optimisation in a regression analysis (Bennion and Bachu, 2008).
for the capillary pressure, relative permeability is either measured or described by semi-
empirical formulations as a function of the wetting-phase saturation. For saturations below
the residual saturations, kr,α is zero for this phase, which means the phase is immobile.
The Brooks & Corey model introduced to describe capillary pressure on the macroscale
(Section 2.2.4) can be used to derive a relative permeability-saturation relation (Burdine,
1953) as:
2+3λBC
λBC
kr,w = Se ,
2+λBC
(2.26)
2 λBC
kr,n = (1 − Se ) · 1 − Se .
Sw −Sw,r
In this approach, the effective saturation used is Se = 1−Sw,r −Sn,r
.
As for the capillary pressure relations, little information is available on measured relative
permeability relations for CO2 -brine systems. Figure 2.12 (right) shows measured relations
for the formations in Canada already mentioned together with a Brooks and Corey model
with parameter settings as discussed before. The relative permeabilities for CO2 have been
extrapolated to a value of one for water saturations lower than the residual water saturations.
This makes the relations useful for simulations involving evaporation of water into the CO2
phase, in which case the water saturation could be lower than the residual value.
Hysteresis
In general, hysteresis is a property of a system that exhibits path-dependence to arrive at
a state of the system. With respect to this study, hysteresis is observed in the capillary
pressure- and relative permeability-saturation relations, which means that the actual value
of pc and kr is dependent on the history of the displacement processes that occurred in the
system. Several reasons are considered to be responsible for this effect (on the microscale),
i.e. a difference in the contact angle Θ for drainage (displacement of the wetting-fluid by
2 Conceptual, Mathematical and Numerical Model 35
the non-wetting fluid) and imbibition (displacement of the non-wetting-fluid by the wetting
fluid) processes (cf. Section 2.2.4), varying pore-throat diameters relevant for the drainage
and imbibition process, and the already discussed effect of entrapped, non-wetting phase
fluid bubbles (cf. Section 2.2.4) causing a change in the residual saturation. In recent years,
various authors (Flett et al. (2005), Spiteri et al. (2005), Juanes et al. (2006), Doughty
(2007), Leicht (2007), Spiteri et al. (2008)) studied hysteretic behaviour by experiments or
numerical investigations and described its effect on various questions of interest related to
CO2 storage in geological formations. However, hysteresis is not included in this study since
literature data on the underlying capillary pressure- and relative permeability-saturation
relations is still sparse.
Thermal Conductivity
where λpm (Sw = 0) is the thermal conductivity of the dry porous medium and λpm (Sw = 1)
is the thermal conductivity in fully water-saturated conditions.
is valid if the chemical potential of a component is identical in all phases. The implica-
tions of this assumption on the mathematical and numerical model will be discussed in
Sections A.1.3 and A.2.4. Local thermodynamic equilibrium can be justified by assuming
slow flow velocities.
2.3.1 Advection
The movement of a fluid phase due to a piezometric head gradient is called advection. The
velocity of a phase in a porous medium can be described by Darcy’s law (Darcy (1856), Bear
(1972)):
%g
v = −k · ∇h, (2.28)
µ
where v is the Darcy velocity, and h is the piezometric head as h = %pg + z. Darcy’s law is
valid for laminar flow indicated by Reynolds numbers smaller than one. The dimensionless
Reynolds number relates inertial to viscous forces and can be calculated for a porous medium
by
vd
Re = , (2.29)
ν
where v is the fluid flow velocity and d is the mean pore diameter.
Darcy conducted his experiments for fully water-saturated conditions. To describe the phase
velocity in a multi-phase system, Darcy’s law is extended to (Helmig, 1997)
vα = −k λα ∇pα − %α g ∇z , (2.30)
where λα is phase mobility, pα is phase pressure, and z is elevation. The mobility is calculated
as
kr,α
λα = . (2.31)
µα
The Darcy velocity describes the flow velocity on the macroscale; thus, it neglects the pore
geometry. To obtain the particle velocity, the Darcy velocity is divided by porosity.
2 Conceptual, Mathematical and Numerical Model 37
qC C
α,a = %α Xα v A. (2.32)
2.3.2 Buoyancy
Fluid flow due to buoyancy forces is already described by the Darcy velocity. It is caused by
density differences within one phase or between different phases. It is an important process
in CO2 storage in geological formations, since the CO2 -rich phase typically has a much lower
density than the water-rich phase (cf. Section 2.2.2), causing the CO2 -rich phase to migrate
upwards (against gravity). To prevent CO2 from rising further, appropriate trapping mech-
anisms are required. This is discussed in Section 1.1.
qh = −λi ∇T A, (2.38)
In a multi-phase system, the temperatures of each phase (including the solid phase) may be
different; hence, heat conduction needs to be described separately for each phase. However,
due to small flow velocities in the problem of interest here which allows enough time for
temperature equilibration among the phases, a single temperature value can be assumed for
all phases. Consequently, a single thermal-conductivity value λpm is used in this study, as
defined in Section 2.2.4.
dissolution
thermal energy
Figure 2.13: Non-isothermal two-phase two-component model concept for the CO2 -water system.
The water-mass fraction in the CO2 -rich phase is however is considered constant in this
concept since it is at least one order of magnitude smaller than the CO2 -mass fraction in the
water-rich phase in the conditions of interest here. Moreover, only few approaches exist in
the literature which describe the change of the CO2 -rich phase properties due to the evapo-
ration/dissolution of water (Bielinski, 2006).
The mass transfer of the CO2 component between the phases is of importance here. Carbon
dioxide dissolves in water, producing a weak acid according to the chemical reaction
where in common notation H2 CO3 is carbonic acid, HCO− 3 is hydrogen carbonate, and CO3
2−
is carbonate. The change of the pH-value and chemical reactions with the rock (porous
medium) or other chemicals (e.g. impurities of the injected CO2 ) are not considered here. Ki-
netic effects are also neglected and the resulting mass fraction of CO2 in the water-rich phase
in equilibrium conditions is dependent on the pressure, temperature, and salt content. Bielin-
ski (2006) compared various approaches by different authors describing an EOS to quantify
the amount of dissolved CO2 . Accordingly, the proposed approach by Duan and Sun (2003)
is used in the following. Figure 2.14 shows the resulting variation of the CO2 mass fraction
in brine with pressure and salinity at different temperatures.
The energy transfer between the phases is included in the model concept by calculating the
water-rich phase enthalpy depending on the CO2 mass fraction (cf. Section 2.2.2).
40 2.5 The Simulation Platform MUFTE-UG
0.05 0.065
mass fraction CO2 in brine [kg/kg]
Figure 2.14: Variation of the CO2 mass fraction in brine with pressure at a constant salinity of
0.1 kg/kg (left) and with salinity at a constant pressure of 100 bar (Bielinski, 2006).
The model follows an Eulerian approach where a set of balance equations is derived for a
fixed control volume. In order to solve this set of balance equations, a number of closure
relations are defined. The set of balance equations is solved numerically for the independent
unknowns, the so-called “primary variables”. Since two levels of complexity are considered
to describe the relevant processes, consequently two types of mathematical and numerical
models are derived in this study. These types of models are called “modules” in the following.
The first module, called 2p-module, considers (only) isothermal multi-phase processes in
porous media (cf. Section A.1). The second module, called 2p2cni-module, considers non-
isothermal multi-phase multi-component processes in porous media (cf. Section A.2). The
sections describing the initial and boundary conditions (Section A.3) and the linearisation
and solution (Section A.5) apply to both modules. The section describing the discretisation
of the differential equations in space and time (Section A.4) is shown exemplarily for the
2p2cni-module. The 2p- and 2p2cni-modules are implemented in the simulation platform
MUFTE-UG, described in the follwing Section.
model concepts and discretisation methods for isothermal and non-isothermal multi-phase
multi-component flow and transport processes in porous and fractured-porous media (Helmig
(1997), Helmig et al. (1998), Class et al. (2002), Assteerawatt et al. (2005)). UG is the ab-
breviation for Unstructured Grid. This toolbox provides the data structures and fast solvers
for the discretisation of partial differential equations based on parallel, adaptive multigrid
methods (Bastian et al., 1997). MUFTE-UG’s special advantages also include the data
structures for unstructured grids, functional parallelisation, especially designed for MIMD
(Multiple Instruction stream, Multiple Data stream) parallel computers and adaptive local
grid refinement.
Institute for Hydraulic Engineering (IWS) Interdisciplinary Center for Scientific Computing (IWR)
The simulation platform MUFTE-UG, especially the 2p- and 2p2cni-modules referred to
here (Sections A.1 and A.2), has been extensively tested in code intercomparison studies
and provides results in good agreement with other commercial and non-commercial codes
(Pruess et al. (2003)). Detailed information on the capabilities of MUFTE-UG to simulate
CO2 storage in geological formations is given in Bielinski (2006).
3 Properties of Potential Geological
Formations∗
To investigate CO2 injection processes in geological formations it is of importance to have a
good knowledge of the range and distribution of the relevant reservoir parameters. These pa-
rameter ranges and distributions form the basis for setting up e.g. typical reservoirs allowing
the definition of a number of random reservoir parameter setups respecting the statistical
characteristics of parameter distributions. Statistical characteristics of reservoir parameters
relevant for CO2 storage on the reservoir scale can be calculated from the U.S. National
Petroleum Council public database (NPC, 1984). The database was developed for the as-
sessment of the United States’ enhanced oil recovery potential in 1984. Today, it is part
of the very comprehensive TORIS (Total Oil Recovery Information System) database. In
this database, 2540 oil reservoirs are listed, accounting for over 64 % of the original oil-in-
place estimated to exist in discovered crude oil reservoirs in the U.S.A.. The entire TORIS
database is unfortunately not available for analysis. Nevertheless, the available public part
comprises about half the total number of reservoirs listed in TORIS.
43
44 3.2 Test on Hypothesised Statistical Distributions
saline formations do not differ fundamentally from those shown in the following for oil and
gas reservoirs derived from the NPC database. Consequently the statistical characteristics
can be used to evaluate CO2 injection and storage in aquifers. The NPC database gives
average parameter values for the entire reservoir. The histograms of porosity, logarithm of
the intrinsic permeability, geothermal gradient, and mean reservoir depth below the surface
are shown in Figure 3.1.
relative relative
frequency frequency
0,20 0,20
0,15 0,15
0,10 0,10
0,05 0,05
0,00 0,00
7,0 17,0 27,0 37,0 47,0 57,0 -16,5 -15,5 -14,5 -13,5 -12,5 -11,5 -10,5
Porosity [%] Log (k/1m²)
relative relative
frequency frequency
0,20 0,20
0,15 0,15
0,10 0,10
0,05 0,05
0,00 0,00
0,01 0,02 0,03 0,04 0,05 0,06 0,07 0,08 0,09 0,10 250 1500 2750 4000 5250
Geothermal Gradient [°C/m] Depth below Surface [m]
Figure 3.1: Histograms data show relative frequency of porosity [%] (top-left), absolute perme-
◦
ability [m2 ] (top-right), geothermal gradient [ mC ] (bottom-left), and reservoir depth
below surface [m] (bottom-right) derived from the NPC database. Lines indicate
normal or log-normal distributions having the same statistical characteristics as the
respective histogram data sets.
Statistical characteristics are given in Table 3.1. Due to data errors in the records, only
a portion of the reservoir data are used in the analysis. A typical error are negative data
entries, e.g. negative porosity. The number of data values used (n) is given for each parameter
in Table 3.1.
Table 3.1: Statistical characteristics for reservoir parameters in the NPC-database: logarithm of
total absolute permeability (Log (k/1m2 )), geothermal gradient (G.grad), porosity (φ),
depth below surface, salinity and dip angle. Stated are the number of values in data
sets (n), minimum (Min), maximum (Max), arithmetic mean (A. Mean), median, 5th
and 95th percentile of the data. Maximum values for absolute permeability, geothermal
gradient, and porosity seem to be unrealistic. However, in the further course of the
study only median, 5th , and 95th percentile values are used (cf. definition of typical
reservoirs in Section 3.4).
between the cumulative probability of the data and the cumulative probability of the hypoth-
esised distribution (e.g. a normal distribution). The test statistics are calculated according
to Equation 3.1.
where Fe(x) denotes the cumulative probability of the data and F (x) the cumulative proba-
√
bility of the hypothesised distribution. The test statistics (multiplied by n) are compared
to a threshold. The threshold depends on a level of significance to be specified. The tested
hypothesised distributions are a normal distribution and a log-normal distribution for all
parameter data sets. See Figure 3.2 for cumulative probability plots of the permeability and
geothermal gradient data sets, the hypothesised distributions and the test statistic (a).
The (null) hypothesis, i.e. that the data sets investigated follow either a normal or log-normal
distribution, was rejected for all reasonable thresholds. Consequently, standard probability
distributions cannot be used further in this study.
1 1
0.8 0.8
cumulative probability [-]
0.4 0.4
~ ~
F(x) F(x)
F(x) F(x)
0.2 a 0.2 a
0 0
-16.0 -15.0 -14.0 -13.0 -12.0 -11.0 0.009 0.059 0.109 0.159
2
Log10 (absolute permeability / m ) Geothermal Gradient [°C/m]
Figure 3.2: Kolmogorov-Smirnov test on absolute permeability and geothermal gradient. The
cumulative probability of a log-normal distribution is shown, having absolute perme-
ability parameters (F (x), left) and a normal distribution having geothermal gradient
parameters (F (x), right), the cumulative probability of the data (Fe(x)) and the dif-
ference between those two (a). Inset: Since the data are a piecewise constant function
and the hypothesised distribution is continuous, a difference can be calculated to the
left (a1 ) and to the right (a2 ) of the considered position x. The larger one is selected
as a.
with depth, and variation of porosity with depth. The correlation coefficients are shown in
the boxes.
Due to the low correlation coefficients, mutual interrelation of investigated parameters can
be rejected. The hypothesized interrelation between absolute permeability and porosity will
be further investigated in Chapter 7.
Case M The median reservoir case assumes median parameters for absolute permeabil-
ity (123 mD† ), geothermal gradient (0.03 ◦ C/m), porosity (20 %), depth below surface
(1524 m) and salinity (0.048 kg/kg). A Brooks and Corey model (Brooks and Corey,
1964) capillary pressure-saturation relation and relative permeability-saturation rela-
tion is assumed employing the input parameter set as given in Table 2.2.
Cases W and C The warm (W) and cold (C) reservoir cases have the same median param-
eters, except for the geothermal gradient; here the 5th and 95th percentile values are
†
1D=9.86932 · 10−13 m2 ≈ 10−12 m2
3 Properties of Potential Geological Formations 47
Figure 3.3: Variation of absolute permeability with porosity, variation of absolute permeability
with depth, and variation of porosity with depth. Interpolated functions of the re-
spective dataset are shown as red lines and correlation coefficients R are given in the
boxes. Correlation coefficients are all rather low, however, an interrelation between
permeability and porosity can be hypothesized.
used (0.018 ◦ C/m and 0.062 ◦ C/m, respectively). This means that these reservoirs are
at the same depth as the median reservoir, but have a higher or lower temperature
compared to the median reservoir.
Cases S and D The shallow (S) and deep (D) reservoir cases have the same median param-
eters, except for depth; here the 5th and 95th percentile values are used (386 m and
3495 m, respectively). The shallow case reservoir is located at a depth of only 386 m
and thus at sub-critical conditions for CO2 . This case might be unrealistic for a large
scale storage attempt. However, I do not want to manipulate the statistical character-
istic of the parameters and moreover consider this case to be important to show the
sensitivity of the results.
Cases V, E, and B The median reservoir parameter setup is combined with measured rel-
ative permeability relations of CO2 -brine systems in the Alberta basin in Canada
(Viking (V)-, Ellerslie (E)-, and Basal (B) Formations) as introduced in Section 2.2.4
and shown in Figure 2.12.
Case P The median reservoir parameter setup is combined with a three times higher capil-
lary entry pressure (pd =30 kPa) in the Brooks and Corey models (Brooks and Corey,
48 3.5 Summary and Conclusion
Case Q The Median reservoir parameter setup is combined with a half the CO2 injection
rate compared to the other cases (the actual value depends on the model setup (1-D
or 3-D) and is given in the relevant sections).
Case K Finally, a combination of the median reservoir parameter setup and an absolute
permeability reduced by one order of magnitude (12.3 mD) yields Case K.
Median Reservoir M 1524 15.47 55.13 0.30 0.05 1025.5 660.7 0.084
Warm Reservoir W 1524 14.95 104.49 0.30 0.05 994.5 316.8 0.078
Cold Reservoir C 1524 15.51 37.43 0.30 0.05 1031.7 805.5 0.088
Shallow Reservoir S 386 4.06 21.58 0.30 0.05 1032.9 98.3 0.015
Deep Reservoir D 3495 35.31 115.00 0.30 0.05 995.2 666.1 0.174
Median Reservoir+V.-kr V 1524 15.47 55.13 0.56 0.044 1025.5 660.7 0.084
Median Reservoir+E.-kr E 1524 15.47 55.13 0.66 0.034 1025.5 660.7 0.084
Median Reservoir+B.-kr B 1524 15.47 55.13 0.29 0.035 1025.5 660.7 0.084
3 Properties of Potential Geological Formations
Median Reservoir+3pd P 1524 15.47 55.13 0.30 0.05 1025.5 660.7 0.084
Median Reservoir+ 2q Q 1524 15.47 55.13 0.30 0.05 1025.5 660.7 0.084
k
Median Reservoir+ 10 K 1524 15.47 55.13 0.30 0.05 1025.5 660.7 0.084
Table 3.2: Parameter settings for typical reservoirs. The three rightmost columns show calculated parameters based on pressure, tem-
perature, and salinity. Cases V, E and B, use median reservoir properties in addition to measured relative permeability data
for the Viking (V), Ellerslie (E) and Basal (B) formations from Alberta, Canada (Bennion and Bachu, 2008). For the case
with index P a tripled entry pressure for the Brooks and Corey capillary pressure-saturation model is used. In Case Q, CO2
is injected at half rate compared to the others. In Case K, the permeability is reduced by one order of magnitude.
49
4 Dimensional Analysis∗
To identify and assess the dominant forces and processes relevant in CO2 storage in geological
formations, a dimensional analysis of the governing equations is performed. Characteristic
values for length, time, velocity, and pressure are introduced in the balance equations in frac-
tional flow formulation. These characteristic values represent typical length scales, typical
timespans, typical flow velocities, and typical pressures which can be found when analysing
the processes of interest, e.g. the CO2 plume evolution. Furthermore, dimensionless num-
bers are defined, representing relations of forces in the system, i.e. capillary, viscous, and
buoyancy forces. The resulting set of equations consists only of dimensionless gradients,
dimensionless numbers and dimensionless functions which depend on phase saturations, rel-
ative permeability relations and fluid property relations. Thus, it is possible to investigate
the individual dimensionless terms of the equations independently of each other by selecting
the characteristic values or dependent on each other by calculating the characteristic values
in numerical simulation experiments. As a basis for investigating the dimensionless terms,
the typical reservoirs are considered as defined in Chapter 3. Two types of reservoirs are
considered, i.e. a 1-D gravitation-free reservoir and a radially symmetric 3-D reservoir. The
aim is to obtain a better qualitative understanding of the influence of the relations of forces
and the selection of the characteristic values on the CO2 plume evolution behaviour, the
resulting storage capacity, and risk.
∗
This Chapter is published in Kopp, A., Class, H. and Helmig, R., Investigations on CO2 storage
capacity in saline aquifers - Part 1: Dimensional analysis of flow processes and reservoir characteristics, Int.
J. Greenhouse Gas Control, 3(3), 263–276, DOI:10.1016/j.ijggc.2008.10.002, 2009.
51
52 4.1 Derivation of Dimensionless Formulation
phase α is defined as
λα
fα = , (4.1)
λ
where the total mobility λ can be calculated according to
X
λ= λα . (4.2)
α
Summing up Equation A.4 over the phases and defining the total velocity vtot as the sum
of the phase velocities (Equation 2.30), the pressure equation can be stated as (Chen and
Ewing, 1997)
∂φ X 1 ∂%α
∇· vtot = − − φ Sα + vα ∇%α − %α qα , (4.4)
∂t %α ∂t
α X
vtot = − λ k ∇pα + fβ ∇(pβ − pα ) − g ∇z fα % α , (4.5)
α
with phases α, β ∈ (w,CO2) and α 6= β. Resolving the phase velocities (Equation 2.30) for
k∇pα , then equating and solving again for the phase velocity yields
The saturation equation is obtained by inserting Equation 4.6 into Equation A.4:
∂Sα
φ + ∇ · fα vtot + ∇ · fα λβ (%α − %β ) k g ∇z + ∇ · fα λβ k ∇(pβ − pα )
∂t
(4.7)
∂φ 1 ∂%α
= − Sα − φ Sα + vα · ∇%α + qα .
∂t %α ∂t
Equations 4.4, 4.5, and 4.7 are the fractional flow formulation for two-phase flow in porous
media. Some assumptions can be made to simplify these equations for further analysis. The
assumption of a rigid rock matrix leads to ∂φ ∂t
= 0. Both fluid phases are assumed to be
∂%α
incompressible, leading to ∂t = 0 and ∇%α = 0. This assumption is not valid for the CO2 -
rich phase in reservoirs at shallow depth, however it is necessary at this point to simplify
the equations. Sources or sinks are not considered; consequently, the right-hand side terms
in Equations 4.4 and 4.7 vanish. In the following, the dynamic viscosity is also assumed to
be constant.
4 Dimensional Analysis 53
Finally, using pw as the independent pressure variable, one is left with a simplified set of
pressure equations
∇· vtot = 0, (4.8)
X
vtot = − λ k ∇pw + fCO2 ∇pc − g ∇z fα %α , (4.9)
α
and one saturation equation each for the water and CO2 phases
∂Sw
φ + ∇ · fw vtot (4.10)
∂t
+∇ · fw λCO2 (%w − %CO2 ) k g ∇z
+∇ · fw λCO2 k ∇pc = 0,
∂SCO2
φ + ∇ · fCO2 vtot (4.11)
∂t
+∇ · fCO2 λw (%CO2 − %w ) k g ∇z
−∇ · fCO2 λw k ∇pc = 0.
The first term in Equations 4.10 and 4.11 is the accumulation term, the second term has
(viscous) advective character, the third term has (gravitational) advective character, while
the fourth term is of (capillary) diffusive nature.
The characteristic length (lcr ) is used to non-dimensionalise the elevation head (Equa-
tion 4.12) and the spatial derivatives (Equation 4.13). When applying lcr to non-dimensionalise
the spatial derivatives, it relates to a length over which a characteristic gradient occurs.
Therefore, this length is alternatively chosen as the system length or the length (width) of
a saturation front. In the first case, the characteristic gradient indicates a saturation or
pressure drop over the system length; in the latter case, it refers to a saturation or pressure
drop over the front length (width). The characteristic time (tcr ) is chosen as the time over
which a characteristic saturation change occurs, and is thus related to the front propagation
velocity. This means that the characteristic total velocity (vcr ) is coupled to the characteris-
tic time and the characteristic length by vcr = φtcrlcr . Hence, one of these characteristic values
tcr , lcr , vcr is always dependent of the other two. The characteristic pressure pcr is used to
non-dimensionalise phase pressures and capillary pressure (Equations 4.15 and 4.16). This
procedure is in agreement with the definition of capillary pressure (Equation 2.23). Hence,
the characteristic pressure is chosen as a measure related to either capillary pressure, e.g. a
capillary pressure drop over the system length or the front length (width), or to phase pres-
sure, e.g. a phase pressure drop over the system length or the front length (width). The
definitions of characteristic values are summarised in Table 4.1.
The choice of characteristic values is further discussed and illustrated in Sections 4.3 and 4.4.
The dimensionless Capillary and Gravitational Numbers relate forces acting on the system,
where Ca relates forces (M·L)/T2 and Gr relates forces per unit length M/T2 (M is mass, L
is length, and T is time). Note that absolute permeability in the dimensionless numbers is
a scalar, which means that permeability is assumed to be isotropic.
A formal relation of capillary forces to gravitational forces can be defined, when dividing the
Capillary Number Ca by the Gravitational Number Gr, this is the Bond number (Bo).
However, the Bond Number is not found in the balance equations (as will be shown in the
following section), therefore it is not discussed any further here.
∇·
b v̂tot = 0, (4.21)
v̂tot = − Ca A ∇p̂
b w − Ca kr,CO2 ∇p̂
b c + Gr A B ∇ẑ,
b (4.22)
where
µCO2
A= kr,w + kr,CO2 , (4.23)
µw
%w fCO2 %CO2 fw
B = (fw − fCO2 ) + + . (4.24)
%w − %CO2 %w − %CO2
∂Sw b
+ ∇ · fw v̂tot + ∇
b · Ca C ∇p̂
b c +∇b · Gr C ∇ẑ
b = 0, (4.25)
∂ t̂
where
C = kr,CO2 fw . (4.26)
P
Using the closure relation α Sα = 1, the saturation equation for the CO2 phase can simply
be derived out of Equation 4.25 as
∂SCO2
−∇
b · (. . .) − ∇
b · (. . .) − ∇
b · (. . .) = 0. (4.27)
∂ t̂
The detailed derivation of Equations 4.21, 4.22, and 4.25 is given in Appendix B.
56 4.3 Analytical Investigations
Storage capacity is, at this point, simply defined as the volume fraction of the reservoir
available for CO2 storage. Hence, it is a dimensionless quantity. Theoretically it ranges
between zero (no storage is possible) to the volume fraction of the entire pore space of the
reservoir, that is average porosity.
10 6 10 6 lcr=20m
M Median M Median
W Warm W Warm
4 4
10 S WM C Cold 10 S WMDP C Cold
S Shallow S Shallow
D Deep D Deep
K Median & k/10 K Median & k/10
10 2 10 2 P Median & 3pd
gravitational capillary
forces forces
Ca [-]
Gr [-]
increase increase
10 0 10 0
viscous viscous
forces D C K forces CK
increase increase
-2 -2
10 10
10 -4 10 -4
10 -6 -10 10 -6 -10
10 10 -9 10 -8 10 -7 10 -6 10 -5 10 -4 10 -3 10 10 -9 10 -8 10 -7 10 -6 10 -5 10 -4 10 -3
vcr [m/s] vcr [m/s]
Figure 4.1: Variation of Gravitational Number Gr versus the characteristic velocity (left) and
variation of the Capillary Number Ca versus the characteristic velocity (right). Letter
coding according to Table 3.2.
As apparent from the definition of Gr, gravitational forces increase with increasing difference
in fluid densities and increasing permeability. Viscous forces increase with increasing char-
acteristic velocity or increasing dynamic viscosity. For the median case reservoir, a reference
equilibrium state between gravitational and viscous forces, as indicated by Gr = 1, is given
for vcr = 8.5 · 10−6 m/s. If Gr becomes significantly larger than one, gravitational forces tend
to dominate over viscous forces and vice versa. The number Gr is very sensitive to variations
in reservoir depth. When comparing Gr calculated from median and from shallow reservoir
properties, an increase in Gr by one order of magnitude for a presumably constant charac-
teristic velocity is observed (the assumption of constant characteristic velocity is discussed
in the following paragraph). This means that gravitational forces increase by a factor of ten
in relation to viscous forces. For the Deep reservoir, fluid property changes are smaller and
Gr decreases by a factor of 2.3 with respect to the Median case value. A similar behaviour
is observed for a warmer / colder reservoir, a large increase in Gr for a warmer reservoir and
a small decrease for a cooler reservoir. This is related to the non-linear behaviour of CO2
density depending on pressure and temperature. With increasing depth (increasing pressure
and temperature) changes in density become smaller (Bachu, 2003). A change in absolute
permeability has a strong influence on Gr. As permeability may vary by several orders of
58 4.3 Analytical Investigations
magnitude, this implies a variation in Gr by the same magnitude. However, with a different
permeability, characteristic velocity will change simultaneously since both are closely coupled
(Equation 2.30), or, depending on the boundary conditions, the pressure gradient will change.
Figure 4.1 (right) shows the variation of Ca versus vcr . It is obvious from the definition of
Ca that capillary forces increase with increasing characteristic capillary pressure or perme-
ability. Viscous forces increase with increasing characteristic velocity, dynamic viscosity, or
characteristic length. When varying vcr , an assumption has to be made for the characteristic
length lcr in order to calculate Ca. Here lcr is taken to be 20 metres which corresponds to
the intermediate front width observed in the simulations (see Section 4.4). For the Median
case, the reference equilibrium state between the forces, indicated by Ca = 1, is given for
vcr = 1.18 · 10−6 m/s. This is a significantly lower velocity than for the reference equilibrium
state between gravitational and viscous forces (Gr = 1). This means that for the Median
case, where capillary forces are at (reference) equilibrium state with viscous forces, the grav-
itational forces dominate over the viscous forces, indicated by Gr = 7.2. This is further
analysed in the discussion of Figure 4.2. The Ca number is less sensitive to variations in
reservoir temperature and depth than Gr. The ratio of capillary forces to viscous forces
increases by a factor of two for the Warm reservoir and by factor of 3.6 for the Shallow
reservoir compared with the Median reservoir at a given characteristic velocity. The same
strong influence of permeability on Gr is found on Ca. The influence of pcr on Ca is strong.
As previously discussed, the characteristic pressure pcr is used to non-dimensionalise capil-
lary pressure and phase pressures. Very little information is available in the literature on
capillary pressure-saturation relations relations for CO2 -brine systems (cf. Section 2.2.4).
Capillary pressures vary over a range of four orders of magnitude for water saturations rang-
ing from zero to one. In CO2 storage scenarios, however, capillary pressures may vary only
by two orders of magnitude, since very low water saturations accompanied by very high cap-
illary pressures occur only in the close vicinity of the injection well. When selecting a phase
pressure as the characteristic pressure, the variation in pressure between a deep and a shal-
low reservoir is approximately one order of magnitude. Therefore, despite a strong influence
of pcr on Ca, the expected changes are not as large as these induced by variable permeability.
M Median&lc=20m
M2 Median&lc=2m
W Warm M S W
C Cold
lc=20m
105 S Shallow
D Deep
K Median & k/10
P Median & 3pd
103
gravitational P K D C M2
101
Gr [-]
forces
increase
viscous
forces
10-1 increase
decreasing vcr
10-3
viscous capillary
forces forces
increase increase
10-5 -5
10 10-3 10-1 101 103 105
Ca [-]
Figure 4.2: Gravitational Number Gr versus Capillary Number Ca for a varying characteristic
velocity ranging from 10−10 m/s to 10−3 m/s. Letter coding according to Table 3.2.
turn, dominate over gravitational forces (lower right sector). In this case gravitational forces
are rather weak, leading to a low buoyancy-induced upward CO2 flow. Viscous forces ideally
should be high, resulting in a cylindrically shaped plume. The stronger the capillary forces,
the more diffuse the displacement front is. Values of both Gr and Ca around unity indicate a
situation in which none of the three forces dominates. Comparing, for example, the Shallow
reservoir to the Median reservoir, the effect of an increase in Gr and an increase in Ca results
in a shift to the top/right. Reducing the permeability of the Median reservoir by a factor
of 10 (Case Median & k/10), results in a shift in the bottom left direction, i.e. in direction
of increasing viscous forces. A reduction of pcr results in a parallel shift to the left, i.e.
decreasing Ca.
Discussing changes in Ca or Gr based on changes in fluid properties or reservoir proper-
ties like permeability, and at the same time assuming a constant characteristic velocity,
might lead to a mis-interpretation. As CO2 density varies strongly between, for example,
a super-critical (Median reservoir) and a sub-critical (Shallow reservoir) case, the charac-
teristic velocity presumably changes as well. Here, a CO2 density decrease by a factor of
6.7 leads to an increase in characteristic velocity of the same factor when applying a simple
volume balance and assuming a constant cross-sectional area. Due to the decrease in CO2
Shallow
viscosity, viscous forces increase only moderately. Consequently, the actual ratio Gr GrMedian
is
only 1.2 (versus 10 when not considering a corresponding increase in characteristic velocity).
Similarly, when reducing the permeability by a factor of 10 and comparing Gravitational
Numbers, the influence of a changing characteristic velocity has to be accounted for.
60 4.3 Analytical Investigations
101 10-1
10-2
C
B
0
10-3
10
10-4
A, B
C
A
10-5
10-1
10-6
10-7
10-2 10-8
0 0.2 0.4 0.6 0.8 1
S CO2 [-]
Figure 4.3: Functions A, B, and C for a relative permeability-saturation relation model after
Brooks and Corey (1964) (employing input parameter setting of λ = 2.0, pd = 10 kPa,
Sw,r = 0.3 and SCO2,r = 0.05) and fluid property relations for the Median reservoir.
The magnitude of these functions, dependent on CO2 saturation, can now be discussed. Note
that the relative permeability endpoints (i.e. relative permeability at residual saturations)
equal one in the Brooks and Corey model assumed here. A CO2 saturation close to zero
refers to the very initial phase of injection or to locations in a (homogeneous) reservoir far
away from the injection well. A high CO2 saturation relates to regions close to the injection
well at later times or to a region below the caprock, where CO2 presumably accumulates.
For CO2 saturations lower than SCO2,r , function A equals the dynamic viscosity ratio of the
two fluids ( µµCO2
w
). With increasing CO2 saturation, A decreases to a minimum and increases
thereafter up to unity at the residual water saturation. Function B is always greater than one
for %CO2 > %2w . It ranges between −1 + %w −%
%w
CO2
at residual water saturation and 1 + %w%−%
CO2
CO2
for residual CO2 saturation (as can be derived from Equation 4.24). Function C varies by
4 Dimensional Analysis 61
Dimensionless gradients should be ideally scaled to unity to allow for estimates of dominant
forces. Characteristic values lcr and pcr have to be found to non-dimensionalise pressure and
elevation head gradients. At different locations within a reservoir, ∇pα , ∇pc , and ∇z can vary.
Pressure gradients vary also with time. But at various points in time or locations within the
reservoir, e.g. right at the CO2 front, these gradients can be estimated. In the vicinity of the
injection well, phase pressure gradients are high, especially at the beginning of the injection
process when relative permeabilities are low. At later times, pressure gradients are low. As
capillary pressure is dependent on saturation, significant capillary pressure gradients (∇pc )
occur where significant gradients in saturation are found. This is the case right at the front
of the CO2 plume. In Section 4.4, it will be shown that it is possible to set either ∇p̂
b w or ∇p̂
b c
equal to unity over either the front width or the system length. The gradient in elevation is
obvious. In the horizontal plane, ∇z is zero and in vertical direction, ∇z is one. According
to the NPC database (see Table 3.1), average reservoir dip angles of up to 25 degrees are
possible, which correspond to ∇z = 0.466.
D
0.6
M
S CO2 [-]
K
0.4 P
S
0.2 V
Q B C W
E
0
0 5E-07 1E-06 1.5E-06 2E-06 2.5E-06
x/t [m/s]
Figure 4.4: Carbon dioxide saturations versus similarity variable (distance from injection point
divided by model time) in a horizontal 1-D reservoir after 4 years modeltime for all
cases. Parameter settings and letter coding as shown in Table 3.2.
The fastest plume evolution is observed for the Shallow reservoir because of low CO2 density.
The front evolves with a velocity of vCO2 = 1.05 · 10−5 m/s (not shown in Figure 4.4).
Similarly, the plume evolution in the Warm reservoir occurs faster than in the Median
reservoir for the same reason. In contrast, plume evolution is slower in the Deep and Cold
reservoir cases. The slowest plume evolution is observed for Case Q, where the injection
rate is halved. The tripled entry pressure (Case P) causes a more diffusive front. A lower
absolute permeability (Case K) causes a slightly slower plume evolution. This is due to the
higher pressure gradient ∇pw necessary to maintain the mass injection rate and at the same
time quasi-constant ∇pc , thus the front is less diffusive (more compact). For Case B (Basal
reservoir relative permeability relation), plume evolution is also markedly slower, with at the
same time much higher CO2 saturations. In this case carbon dioxide relative permeability is
significantly lower than for the Median case for CO2 saturations lower than 0.7. Viking and
4 Dimensional Analysis 63
Ellerslie reservoir relative permeability relations produce almost identical results, a faster
plume evolution with lower CO2 saturations.
To calculate Ca, characteristic values vcr , lcr , and pcr have to be defined. The character-
istic velocity vcr is estimated from the analytical solution of the Buckley-Leverett problem
(Buckley and Leverett, 1942), analysing solely the fractional flow functions. This procedure
is depicted in Figure 4.5. Note that capillary effects are neglected, and a hyperbolic problem
is assumed for the determination of vcr .
1-S w,r
fCO2 [-]
S CO2,r
S C02t
0.5
tangent to
fCO2
0
0 0.5 1
S CO2 [-]
Figure 4.5: Construction of a sharp front (modified after LeVeque (1992) and Helmig (1997)) by
taking the Brooks and Corey model fractional flow function as an example. The front
is reported to travel with a speed vt proportional to the derivative of a linear function
starting at (SCO2 = SCO2,r and fCO2 = 0) and being tangential to the fractional flow
t
function fCO2 . The saturation of the tangential point is called SCO2 further on.
The linear correlation coefficient (Pearson r, reflecting the degree of linear relationship be-
tween two variables) is r = 0.9998, obtained from the evolution velocity of the front with a
t
saturation of SCO2 observed in the simulation experiments and from the analytically derived
q
vt times the injection rate q divided by porosity φ and CO2 density %CO2 , i.e. vt φ %CO2 . This
means that there is an almost perfect positive linear relationship between the two variables.
This allows an estimation of the characteristic velocity without using model results. As out-
lined in Section 4.3.3, dimensionless gradients ∇p̂
b w and ∇p̂
b c should ideally be scaled to unity.
This can either be obtained for the front or for the entire system. When considering a front
approach, lcr refers to the front width and pcr refers to the pressure drop over this front width
(water pressure or capillary pressure). The front width lcr is defined as the distance between
t
the farthest point the front has travelled and the point where SCO2 is observed. When con-
sidering a system approach, lcr refers to the system length (in this case the model domain)
and pcr refers to the pressure (water/capillary) drop over the system length. Consequently,
four options arise:
b c =! 1 over the front width. See Figure 4.6 for the Median case construction sketch.
1. ∇p̂
64 4.4 Numerical Investigations
• lcr is defined as the front width, i.e. distance between the point where CO2
t
saturation is zero and the point where CO2 saturation is equal to SCO2 .
• pcr is defined as the capillary pressure drop over this front width:
pcr = pc,SCO2
t − pc,front .
Equations 4.13 and 4.15 yield:
p t −pw,front pw,S t −pw,front
b w = w,SCO2
∇p̂ · plcrcr = CO2
.
lcr pcr
x (S CO2=S CO2t )
0.6
x (S CO2=0.0)
S CO2 [-]
0.4 lcr
M
S CO2t
0.2
0
0 50 100 150 200 250 300
x [m]
18000 x (S CO2=S CO2t )
16000 x (S CO2=0.0)
lcr
pc [Pa]
14000
M t
pc (S CO2 )
12000 pcr
10000 pc,front = pc,entry
0 50 100 150 200 250 300
x [m]
!
Figure 4.6: Selection of characteristic values to obtain ∇p̂
b c = 1 over the front width for the Median
case.
Since phase pressure and capillary pressure are strongly coupled via Equation 2.23, it is
possible to scale gradients ∇p̂
b c and ∇p̂
b w in options 1 and 2 (front approach) simultaneously
to unity.
To calculate v̂tot from Equation 4.22, A, B and kr,CO2 are additionally evaluated at the
t
saturation SCO2 . When comparing the resulting CO2 phase velocity with the numerically
calculated CO2 (Darcy) velocity, the linear correlation coefficient (Pearson r) for the front
approach is r = 0.9871 and r = 0.9779 for the system approach. Note that the results are
identical when either capillary pressure drop or water pressure drop is chosen to be pcr . They
differ depending on the definition of lcr . When looking at the Capillary Numbers calculated
b c =! 1, the order for the cases in Table 3.2 is as follows:
by the front approach with ∇p̂
CaB = 0.80 > CaQ = 0.45 > CaP = 0.38 > CaD = 0.29 > CaC = 0.19 > CaM = 0.18
> CaW = 0.10 > CaV = 0.10 > CaE = 0.09 > CaK = 0.03 > CaS = 0.02
The order of Capillary Numbers closely follows the observed order of evolution velocities
except for the lower permeability case and the tripled capillary pressure. Modification of
permeability k and capillary entry pressure pd alters Ca by orders of magnitude, although
the plume evolution velocities do not change much in the 1-D model approach. The left
boundary condition assumes a constant injection rate, therefore the pressure increase is high
for reduced permeability (k) and diffusive spreading loses influence. In three dimensions, the
pressure increase would be much smaller and rapidly decline with distance from the injection
well. The influence of modified relative permeability, temperature, depth and injection rate
is represented well by Ca.
t
In conclusion, for this 1-D gravitation-free reservoir, SCO2 (the saturation of the tangential
point of the fractional flow function) is used as an estimate for average CO2 saturation. The
evolution velocity of the CO2 front is estimated by dimensional analysis. This is a prerequi-
site for the estimation of storage capacity and risk. The 1-D gravitation-free reservoir shows
that relative permeability plays an important role in plume evolution behaviour. Never-
theless, knowledge of phase pressure increase and the capillary pressure-saturation relation
is necessary, which is difficult to obtain without simulations and further investigations on
reservoir properties.
66 4.4 Numerical Investigations
Aw CO2-front
bw
mCO2
rw rb
rf
r
Figure 4.7: Sketch of a radially symmetric domain. The front velocity at the top of the domain
(at distance rf ) is higher than the velocity at the bottom of the domain (at distance
rb ) due to gravity segregation. Symbols: rw indicates injection well radius, bw width
of injection well sector and Aw the injection well surface.
Due to gravity segregation, the front velocity at the top of the domain (at distance rf ) is higher
than at the bottom (at distance rb ) of the domain. This behaviour cannot be approximated
by a simple volume balance approach. Roughly, the squares of the radial distances are
linearly proportional to time (Barenblatt et al., 1990). For example in the Median case,
the front at the bottom of the domain even stops entirely after 3.7 years. This means that
under these circumstances, regions at the bottom of the domain that are further away than
approximately 150 m from the injection well cannot be accessed for storage. Figure 4.8 shows
CO2 saturation iso-lines for all reservoir setups after 3.85 years of injection and consequently
an injected total CO2 mass of 3.85 Mt CO2 (1.925 Mt CO2 for Case Q).
Here, the focus is on the evaluation of forces acting in the reservoir, and thus on the as-
sessment of dimensionless numbers. To non-dimensionalise the equations and calculate the
dimensionless numbers Ca and Gr, characteristic values vcr , lcr and pcr have to be defined
once again. These values change with time and location in the reservoir. Characteristic
velocity vcr is defined by averaging CO2 Darcy velocities at every timestep in the simula-
tions. Performing mass averaging, the element value (e.g. euclidian norm of the velocity) is
weighted by the CO2 phase mass fraction of the element with respect to the total CO2 phase
4 Dimensional Analysis 67
100
W S
B
Q V
z [m]
P
D
50 M E
C
K
0
0 500 1000 1500 2000
r [m]
Figure 4.8: Carbon dioxide saturation 10% iso-lines developing for relative permeability-
saturation relations given in Figure 2.12 and reservoir cases shown in Table 3.2 after
3.85 years model time. Letter coding according to Table 3.2.
elements
X massCO2,j · || vCO2,j ||e
vcr = Pelements . (4.28)
j=1 i=1 massCO2,i
Two options are examined, ∇p̂ b c =! 1 over the front width and ∇p̂
b c =! 1 over the system length.
Figure 4.9 shows the resulting Capillary Numbers Ca versus Gravitational Numbers Gr. For
each case, the start of injection is marked by a square. Time is reflected by the length of the
curves in Figure 4.9 (however, curve length is not a scale for time). The lengths of the curves
vary, since injection is stopped when the CO2 plume spreads a radial distance of 1 kilometer
from the injection well. For example, the curve for Case K (Median & k/10) describes a
timespan of 7.48 years whereas Case S represents only a timespan of 0.37 years. The Median
case timespan is 3.85 years.
For the estimate, where ∇p̂b c =! 1 over the front width (Figure 4.9, left), Gr is increasing with
increasing time, whereas Ca is rather constant. The Gravitational Number Gr is increasing
with time since vcr is decreasing and all other values in Gr are constant. The Capillary
Number Ca varies only marginally since lcr (selected here as the front width at 32 of the
reservoir height) is increasing with time due to capillary diffusion and vcr is decreasing at
the same rate; these effects cancel each other out. The small oscillations in Ca are due to lcr
being a step function due to numerical discretisation. Note that the characteristic pressure
pcr is constant in time for both options; i.e. capillary pressure drop over front width or system
length. With increasing time, the gravitational forces become more and more important and
finally dominate.
For the estimate where ∇p̂ b c =! 1 over the system length (Figure 4.9, right), both Gr and Ca
are increasing with time. This is due to the fact that vcr is decreasing whereas lcr (selected
as the system length) is constant. The characteristic pressure pcr is slowly increasing here.
68 4.4 Numerical Investigations
Median Median
Warm Warm
Cold Cold
Shallow Shallow
Deep Deep
102 102
Median & Viking Rel.Perm Median & Viking Rel.Perm
Median & Ellerslie Rel.Perm Median & Ellerslie Rel.Perm
Median & Basal Rel.Perm Median & Basal Rel.Perm
Median & 3pd Median & 3pd
Median & q/2 Median & q/2
Median & k/10 Median & k/10
Gr [-]
&
increasing lcr
gravitational gravitational
forces forces
increase increase
100 100
viscous viscous
forces forces
increase viscous capillary increase increasing time viscous capillary
forces forces & forces forces
increase increase decreasing vcr increase increase
10-1 -4 10-1 -4
10 10-3 10-2 10-1 100 101 10 10-3 10-2 10-1 100 101
Ca [-] Ca [-]
Figure 4.9: Variation of the Gravitational Number Gr versus Capillary Number Ca for a radially
symmetric domain derived from simulation experiments. Left: Gr versus Ca where
b c =! 1 over the front width. Right: Gr versus Ca where ∇p̂
∇p̂ b c =! 1 over the system
length.
This means that, with increasing time, and according to this definition of Ca under the
b c =! 1 over the system length, both gravitational and capillary
given assumption that ∇p̂
forces become more and more important and finally dominate over the viscous forces.
When looking at Figure 4.9, it might be confusing that, depending on the choice for char-
acteristic values, the importance of capillary forces varies in time. This behaviour relates
to the scale on which the problem is considered. When looking at the system at a reservoir
scale (lcr = system length), capillary forces are initially weak and viscous forces dominate.
Over time, viscous forces become weaker (due to a reduction of vcr ) and capillary forces
gain importance. Due to this reasoning, reservoir engineers, mostly interested in average
reservoir pressures and production rates, often neglect capillary forces when viscous forces
are strong (field in production). When looking at the system on a process dependent scale
(lcr = front width), capillary pressure is of equal importance during the entire process, since
the process dependent length scale increases at the same rate as the velocity decreases, lead-
ing to constant viscous forces. Considering the system on a process-dependent scale allows
for more detailed and physically motivated interpretations, but requires better knowledge of
the system.
The choice of characteristic values leading to estimates of Ca and Gr as shown in Figure 4.9
has now been physically justified by the processes in the reservoir and the dimensionless
numbers are interdependent now. In the initial analytical investigation (Figure 4.2), this
was not the case because the analysis of Ca and Gr was mathematically motivated and the
characteristic values were independent of each other.
4 Dimensional Analysis 69
Similarly, the effect of the dimensionless numbers on risk can be described. A higher Gr is
expected to result in an increase of risk, since an increase of gravitational forces leads to
enhanced gravity segregation. Hence, a larger fraction of the injected CO2 accumulates in
the upper part of the aquifers, i.e. below the caprock, which is potentially risky. Accordingly,
a lower Gr stands for a lower risk due to an increase of viscous forces, leading to a more
compact plume, and thus larger portions of the injected CO2 are found in deeper parts of the
reservoir. A higher Ca describes an increase in capillary forces in relation to viscous forces.
This results in lower CO2 saturations behind the front due to increased capillary dispersion.
Lower CO2 saturations result in lower risk, since potential leakage rates are lower, if a fault
etc. is encountered. Counteracting on the lower risk estimate is the simultaneously larger
plume extent (at constant density), which increases the risk. However, a stronger influence
of (capillary) diffusive or dispersive processes is expected to result in lower risk since this
generally results in higher fractions of total injected CO2 mass dissolved in brine and thus
leads to less potential leakage of free-phase CO2 .
70 4.5 Summary and Conclusion
• The 1-D simulation study neglects gravity effects and thus reflects the influence of the
Ca number on plume evolution velocities. It was shown that the order (ranking) of
the Ca numbers resembles the order of plume evolution velocities and that the average
CO2 saturation can be estimated by analysing the fractional flow function. This is
a prerequisite when estimating storage capacity. However, particularly the influence
of varying injection rate, capillary pressure and permeability is quantitatively less
significant. The reason for this is that gravity segregation is neglected due to the 1-D
assumption.
• In the 3-D simulation study, the influence of gravity segregation is included and the
development of both the Ca and the Gr numbers was shown. It is possible to order
geological reservoirs qualitatively according to their plume evolution behaviour using
the dimensionless numbers Ca and Gr. The Gravitational Number Gr appears to have
a stronger influence than the Capillary Number Ca. However, it is not clear how the
simultaneous variations of Gr and Ca quantitatively affects storage capacities and risk
estimates. One should expect that Ca numbers further lose influence when the system
dimensions become larger. Nevertheless, a low ratio of gravitational to viscous forces
(low Gr) - and to some extent also a high ratio of capillary to viscous forces (high Ca)-
possibly leads to high CO2 storage capacity and low risk.
71
72 5.1 Discussion of Storage Capacity
pyramids have been proposed, reflecting the multi-faceted aspects of CO2 storage. This
study always refers to the Techno-Economic Resource-Reserve pyramid.
Figure 5.1: Techno-Economic Resource-Reserve pyramid for CO2 storage capacity in geological
media after Bachu et al. (2007).
This pyramid consists of four levels of capacity estimates. The theoretical capacity VCO2,t
[m3 ] represents the first level and assumes that the entire pore volume of a formation minus
the residual water-rich phase saturation is accessible to store CO2 .
ZZZ
VCO2,t = φ(1 − Sw,r ) dxdydz, (5.1)
where Sw,r is residual water-rich phase saturation. The theoretical capacity, which is a
volume, represents the size of the resource pyramid. The second level, the effective capacity,
excludes portions of the theoretical capacity which are not accessible for storage due to
geological and engineering reasons. The third level, the practical capacity, is a subset of the
effective capacity which also considers economic, legal and regulatory, infrastructural, and
general economic aspects. The top level of the pyramid, the matched capacity, is a subset of
the practical capacity. It is obtained by detailed matching of large stationary CO2 sources
with adequate geological storage sites.
To calculate the mass of CO2 that can effectively be stored, Meff , effective capacity is multi-
plied by CO2 density. This can be a delicate task, since pressure and temperature vary during
the CO2 injection phase depending on many reservoir and process controlled parameters (i.e.
absolute and relative permeability, reservoir boundaries, injection strategy, etc.). Note that
the effective capacity, CO2 density, and consequently Meff are generally time dependent.
The aim of this section is to provide more insight into the effective capacity estimates and the
effective CO2 mass that can be stored (Meff ), based on investigations on reservoir properties
of potential formations (Section 3). No reference to a specific reservoir, basin, or region is
5 Analysis of Storage Capacity 73
made. Instead, the general range of effective storage capacity and the mass which can be
stored therein are investigated, together with processes and parameters which may increase
or decrease them. Several attempts have been made in the literature to estimate effective
storage capacity and Meff , which include the models described in the following.
where Ci is the intrinsic capacity coefficient [-], Cg is the geometric capacity coefficient
[-], Ch is the heterogeneity capacity coefficient [-] and φavg is average formation porosity
[-]. Figure 5.2 illustrates the definition of the capacity coefficients. The intrinsic capacity
coefficient Ci has contributions from the fraction of pore space that is occupied by the CO2 -
rich phase (Ci,CO2 ) and the fraction of the pore space that CO2 , dissolved in the water-rich
phase, would occupy if it was converted to the CO2 -rich phase (Ci,w ). The coefficients Ci,CO2
and Ci,w add up to Ci .
Doughty et al. (2001) estimate Ci,CO2 by a Buckley-Leverett type analysis (Buckley and
Leverett, 1942) as
Ci,CO2 ∼
= SCO2 , (5.3)
where SCO2 is the average CO2 -rich phase saturation [-]. The coefficient Ci,w is estimated by
%w
Ci,w ∼
= Sw · XwCO2 · , (5.4)
%CO2
where Sw is the average water-rich phase saturation, XwCO2 the average mass fraction of CO2
dissolved in the water-rich phase. Parameters and variables Sw , XwCO2 , and % are averaged
behind the CO2 front.
The geometric capacity coefficient Cg accounts for partially penetrating wells, gravity segre-
gation and dipping aquifers which reduce the theoretical storage capacity. It is defined here
as the volume fraction of the entire pore space, occupied by the CO2 -rich phase, divided by
the entire available pore space (limited by the reservoir size, for example a spill point, etc.).
The heterogeneity capacity coefficient Ch accounts for heterogeneities in absolute perme-
ability, leading to a further reduction or increase in accessible storage capacity. Note that
74 5.1 Discussion of Storage Capacity
these capacity coefficients cannot always be estimated individually. For example, in a three-
dimensional, heterogeneous reservoir, only the product of Cg and Ch can be estimated. For
simplified problems, for example, in one dimensional gravitation-free reservoirs, Ci,CO2 and
Ci,w can be estimated individually and Cg and Ch are unity by definition.
Figure 5.2: Sketch showing estimation of storage capacity coefficients depending on model com-
plexity (after Doughty et al. (2001)); dark-grey colour indicates water-rich phase and
light-grey colour CO2 -rich phase. Left: Estimation of intrinsic storage capacity coef-
ficient (Ci ) in a gravitation free, homogeneous reservoir (Cg = Ch = 1). Saturations
etc. are averaged within the light-grey area. Centre: Estimation of the intrinsic (Ci )
and geometric (Cg ) storage capacity coefficients in a homogeneous reservoir (Ch = 1).
Coefficient Cg is equivalent to the ratio of the light-grey area to the dark-grey area.
Right: Estimation of the intrinsic (Ci ), geometric (Cg ), and heterogeneity (Ch ) stor-
age capacity coefficients in a heterogeneous reservoir. Note that in the latter case,
only the product of Cg and Ch can be estimated (that is the ratio of the light-grey
area to the dark-grey area), since individual contributions can not be identified.
Doughty et al. (2001) do not give further insight into calculation of the effective CO2 mass
that can be stored in this volume fraction of the reservoir (Meff ).
where VCO2,e is the effective storage volume [m3 ] and Cc is the capacity coefficient [-]. Here,
the porosity φ is already included in the calculation of VCO2,t . The capacity coefficient Cc
ranges between zero and one, i.e. between no storage capacity and usability of the entire
theoretical capacity. Bachu et al. (2007) state that Cc incorporates the cumulative effects of
heterogeneities in absolute permeability, CO2 buoyancy and sweep efficiency.
5 Analysis of Storage Capacity 75
To calculate Meff , Bachu et al. (2007) define two limits for the pressure to calculate an
upper and a lower estimate of the CO2 density. The lower estimate of the CO2 density is
calculated by considering initial reservoir pressure conditions, pi . The upper estimate of the
CO2 density is calculated by considering a maximum pressure, pmax , which is the lowest
of the maximum pressure allowed by regulatory agencies in order to avoid rock fracturing,
and the threshold entry pressure of the caprock. For temperature, the authors define initial
reservoir temperature. These definitions lead to a range of Meff as given in Equation 5.6.
Meff,min = %CO2 (pi , Ti ) · VCO2,t · Cc ≤ Meff ≤ %CO2 (pmax , Ti ) · VCO2,t · Cc = MCO2,max . (5.6)
Bachu et al. (2007) state that currently no values for the site-specific capacity coefficient Cc
are provided in the literature.
In this case, the geometric capacity coefficients Cg,CO2 and Cg,w for the CO2 -rich phase
and the water-rich phase differ. Consequently, CO2 -rich phase saturation in Equation 5.3
(Ci,CO2 ) is averaged behind the CO2 -rich phase front, whereas the properties Sw , XwCO2
and % in Equation 5.4 (Ci,w ) are averaged behind the dissolved CO2 front. Probst (2008)
investigates storage capacity by using this approach for multiple realisations of heterogeneous
permeability fields in a realistic reservoir at shallow depth.
However, in this study a continuous injection into a homogeneous reservoir is considered and
it is not necessary to split Cg and Ch , hence the effective storage capacity coefficients C are
identical using either definition.
To calculate Meff , based on Doughty’s model, the effective capacity coefficient C is multiplied
by CO2 density and the total geometric reservoir volume (bulk volume) Vtotal . It is necessary
76 5.2 Numerical Investigations
CO2 (dissolved)
Water
Figure 5.3: Estimation of storage capacity coefficients according to Equation 5.7 for long-term
investigations, that is when the pore volume containing dissolved CO2 (dissolved
share indicated by orange plus yellow coloured area) is considerably larger than the
pore volume containing CO2 -rich phase (volumetric share indicated by yellow coloured
area).
to multiply by the total geometric reservoir volume because average reservoir porosity and
residual water-rich phase saturation (inherently) are included in the coefficient C. The
effective CO2 mass stored Meff is calculated according to Equation 5.8.
This relationship allows portability and comparison of results obtained by the different mod-
els.
radially symmetric domain. The parameter setups for the various cases are identical to those
previously defined (cf. Table 3.2).
To derive capacity estimates, it is necessary to limit the reservoir volume at some point.
This serves to calculate the total geometric reservoir volume Vtotal (cf. Equation 5.8). The
capacity estimates in the following are given at the point in time when the CO2 phase reaches
a spill point and leakage occurs. This “virtual spill point” is defined at one kilometre distance
from the injection well. The model domain extends beyond the leaky well distance to limit
unwanted influence of boundary conditions on model results.
Table 5.1: Storage capacity coefficients for the 1-D gravitation-free reservoir. The coefficient Ci,w
is calculated according to the solubility model of Duan and Sun (2003), but it is
not included in the numerical simulations and consequently it is not included in the
calculation of C. Coefficients Cg and Ch equal unity by definition here. Porosity φ,
necessary to calculate the effective capacity C is 0.2. Note that C refers to the total
geometric volume of the reservoir. To convert C to pore space volume utilised, one
needs to divide C by porosity. The effective mass Meff [-] indicates CO2 mass injected
when reaching the virtual spill point of the reservoir, normalised by the Median case
mass.
and surprising at first sight. But only a normalised mass of 0.48 can be stored, since the
CO2 density is much lower. Thus, the effective capacity coefficient C as defined here, is not
a sufficient measure when analysing different reservoirs. Instead, the resulting Meff has to
be considered.
The Cold and Deep reservoir cases have a high effective mass (Meff ) of CO2 that can be stored
due to high CO2 density. The Median reservoir with Basal formation relative permeability
relations has a higher Meff than the Median reservoir case, since in this case the relative
permeability to CO2 is much lower and therefore average CO2 -rich phase saturation Ci,CO2
is higher. The Median reservoir with either tripled entry pressure, halved injection rate, or
reduced permeability yields similar Meff as the reference case does.
Lower Meff is achieved in the Warm and Shallow reservoir cases due to low CO2 density. The
high value of Ci,w in the Shallow reservoir case (0.234) is also due to the low CO2 density at
these pressure and temperature conditions. A value of 0.234 for Ci,w means that the amount
of CO2 that could be dissolved in brine would occupy a fraction of 0.234 of the pore space
5 Analysis of Storage Capacity 79
when converted to the CO2 -rich phase. The amount of CO2 that could be stored in the brine
phase (dissolved), is almost identical to the amount that could be stored in the CO2 -rich
phase (0.245).
The Median reservoir with Viking or Ellerslie relative permeability relations has a lower
Meff than the Median reservoir case due to lower average CO2 -rich phase saturation. In
these cases the very high residual water-rich phase saturation together with a slightly higher
relative permeability of CO2 in the important range of SCO2 between zero and 0.3 lead to
lower average CO2 -rich phase saturation. The lower relative permeability of the water-rich
phase seems not to be of importance here.
The 1-D gravitation-free reservoir shows that relative permeability plays an important role
for estimating storage capacity. It influences the capacity estimates to a similar extent
as the entire range of reservoir properties like geothermal gradient and depth. Effects of
density are nicely demonstrated. However, the 1-D approach reveals effects of changing
entry pressure, injection rate and absolute permeability only to some extent. Particularly, a
different absolute permeability is expected to have a huge influence on capacity estimates.
This will be shown in the following section.
0.4 M D
W Q
P
0.3 C
Ci,CO2 [-]
K
E
V
0.2
S
0.1
0
0 1 2 3 4 5 6 7 8
Time [yr]
Figure 5.4: Variation of the intrinsic CO2 -rich phase capacity coefficient Ci,CO2 (∼
= SCO2 ) in time
for various reservoirs. Letter coding according to Table 3.2.
Figure 5.5 shows the geometric capacity coefficient Cg versus time. Due to the virtual spill
point, the reservoir size is fixed and therefore Cg increases until the spill point is reached by
the CO2 plume.
A projection of the curves in Figures 5.4 and 5.5 to the r-axis gives the time since the start
of injection and the non-normalised Meff (in Case M this is 3.85 Mt CO2 ). Case Q is the
only exception here. Since the injection rate is halved, the r-axis values need to be halved
as well to represent Meff .
To summarise, Table 5.2 shows capacity coefficients at the point in time when the virtual
spill point is reached. The point in time when the virtual spill point is reached, varies for the
different cases. Therefore, Table 5.2 requires a closer look to be comparable to Figure 4.8.
Figure 4.8 shows CO2 -rich phase saturation iso-lines after 3.85 years of injection, i.e when the
CO2 plume reaches the virtual spill point in the Median reservoir case. The total geometric
volume of the reservoir up to the virtual spill point is
According to Equation 5.8 and the definition of the C value from Table 5.2, the effective
stored mass Meff for the Median case is
Since 1 Mt CO2 per year is injected and the virtual spill point is reached after 3.85 years, the
effective mass Meff for the Median case must be 3.85 Mt CO2 .
In the Median reservoir case (M), a fraction of 0.279 of the pore space (Cg ) is used for storage.
Multiplying Cg with the average CO2 -rich phase saturation (Ci,CO2 ) of 0.3325 and porosity
5 Analysis of Storage Capacity 81
0.7
0.6
0.5 V
E
0.4
Cg [-]
K
C
P
0.3 B
S
0.2 W Q
0.1 M D
0
0 1 2 3 4 5 6 7 8
Time [yr]
Figure 5.5: Variation in time of the geometric capacity coefficient Cg for various reservoirs. Letter
coding according to Table 3.2.
gives an effective capacity coefficient of 0.01855. To convert the effective capacity coefficient
C (which refers to total reservoir volume) to pore space volume, one needs to divide C by the
porosity, yielding approximately 0.093. Generally, the range of effective capacity coefficients
for the radially symmetric reservoir is larger than for the 1-D gravitation-free reservoir. The
coefficient C ranges in this 3-D study from 0.0117 (Case S) up to 0.036 (Case K). The same
cases also mark the smallest and largest value of Meff ; from 0.09 up to 1.94. This means
that in a reservoir at Shallow case conditions (Case S) less than 5 % of CO2 can be stored
compared to an equally-sized reservoir having a smaller permeability (Case K). Note that
the Shallow case assumes sub-critical conditions for the CO2 phase.
The Median reservoir case with reduced permeability (Case K) shows a very high effective
stored mass (Meff ) due to the high geometric capacity coefficient Cg . In other words, Case
K refers to a more cylindrical plume evolution.
Median reservoir cases with relative permeability relations from the Viking (Case V), Ellerslie
(Case E), or Basal (Case B) formations show all higher Meff than the Median case. The
average CO2 -rich phase saturation (Ci,CO2 ) in each case is comparable to the 1-D experiments
and is presumably caused by the relative permeability of CO2 (for Cases V and E, the average
CO2 -rich phase saturation (Ci,CO2 ) is lower compared to the Median case, for Case B it is
higher). However, the geometric capacity coefficients Cg are higher in all cases (V, E, and B)
compared to the Median case, which is presumably due to the lower relative permeability of
the water-rich phase at high water-rich phase saturations, i.e. Sw > 0.85 (at the front). Since
cases V and E have even lower relative permeability of the water-rich phase in the considered
range than Case B, these cases reach even higher Cg . As a result, when calculating Meff , Case
B stores 188 % of the mass that is stored in the Median reservoir. Cases V and E store 129 %
and 145 % respectively. This shows that very different relative permeability relations (see
82 5.2 Numerical Investigations
Table 5.2: Storage capacity coefficients for a radially symmetric reservoir. The coefficient Ch
equals unity by definition here. Porosity φ = 0.2, necessary to calculate effective
capacity C. The coefficient Ci,w is calculated according to the solubility model of Duan
and Sun (2003), but not included in the numerical simulations and consequently not
included in the calculation of C. The effective Mass Meff [-] indicates CO2 mass injected
when reaching the virtual spill point of the reservoir, normalised to the Median case.
A value Meff of 1.0 equals 3.85 Mt CO2 , i.e. in the Median case this mass is injected
before the virtual spill point is reached. Note that C refers to the total geometric
volume of the reservoir. To convert C to pore space volume utilised, one needs to
divide C by porosity (here 0.2).
Figure 2.12) can lead to high Meff . This shows that not only the residual water- and CO2 -rich
phase saturations are of importance when estimating storage capacity, but also the shape of
both branches of the relative permeability relation (exponent λ in the Brooks-Corey model,
cf. Table 2.2 for values for Cases V, E, and B). Case V is shown in brackets in Table 5.2
since fingering occurred here caused by an unstable front. This matter is not pursued further,
since including dissolution of CO2 in brine will reduce fingering. In a realistic heterogeneous
permeability case, the distribution of permeability rather than hydrodynamic instability
dominates fluid displacement (Garcia and Pruess (2003), Johannsen et al. (2006)).
Warm and Shallow reservoir cases show a capacity coefficient C lower than the Median case
together with a low CO2 density, leading to very low Meff . In contrast, Cold and Deep
reservoir cases show a capacity coefficient C higher than in the Median case along with a
high CO2 density, leading to high Meff . In a radially symmetric domain, the tripled capillary
entry pressure leads to an 11 % increase of effective stored mass. This is due to the increase
5 Analysis of Storage Capacity 83
in the geometric capacity coefficient. The CO2 plume is more spread out in this case. An
effective stored mass of 0.77 for Case Q indicates that a high injection rate is advantageous.
Low viscous forces lead here to stronger gravity segregation and a rather inefficient utilisation
of the available pore volume (Cg ).
The qualitative dependency of storage capacity on dimensionless numbers suggested in Sec-
tion 4.4.3 can now be proven by comparison with calculated capacity coefficients.
Applying the upper estimate for the storage capacity coefficient C in this study of 0.036
(compare to Case K) to the lifetime emissions of a typical coal-fired power plant yields a CO2
plume footprint area of 10.51 km2 . This is a plume with a radius of 1.83 km. Assumptions
are a production of 1 Mt CO2 per year over a 25 year period, CO2 density of 660.7 kg/m3
(i.e. density of the Median reservoir case) and 100 m formation thickness.
84 5.3 Summary and Conclusion
• The investigations show a clear tendency that deep, cold and/or low-permeability
reservoirs are generally more favourable for the efficient utilisation of available stor-
age volume than shallow, warm, and/or high-permeability reservoirs. This conclusion
confirms the qualitative analysis of Bachu (2003) for screening and selecting sites for
CO2 storage. An increase in the injection rate also revealed some potential to increase
storage capacity, since higher viscous forces allow for better spreading of the injected
CO2 . In the 3-D simulations, the difference in storage capacity between the most
favourable and least favourable reservoirs was larger than in the 1-D study. The result
of the 3-D study identifies the Median reservoir with reduced permeability, the Median
Reservoir with Basal formation relative permeability relations and the Cold reservoir
(in this order) as the most favourable reservoirs for CO2 storage. Least favourable for
CO2 storage are the Warm reservoir and the Shallow reservoir.
• In the 3-D simulations, the capacity coefficient C ranges from 0.0117 (Case S) to 0.036
(Case K). The same cases also have the smallest and largest value of Meff (between
0.09 up to 1.94, normalised to the median case). This means a shallow reservoir can
store less than 5 % of the CO2 that can be stored in an equally-sized reservoir with
a 10 times lower permeability at median depth. Note that the Shallow case assumes
sub-critical conditions for the CO2 phase.
relations, both branches kr,CO2 and kr,w , together with the residual water- and CO2 -rich
phase saturations are of high importance. This study indicates that kr,CO2 together
with the residual water-rich phase saturation have a major influence on Ci,CO2 , whereas
kr,w has a major influence on Cg . As already indicated in Section 4.4, the relative
permeability-saturation relations influence the capacity estimates to a similar extent
as the entire range of reservoir properties like geothermal gradient and depth. Since
few data are available on these relations, particularly for CO2 -brine mixtures, this
indicates a great need for further research.
• The results have shown that in all 3-D cases less than 3.6 % of the total geometric
reservoir volume could be used by the injected CO2 before the spill point at one kilo-
metre distance was reached. This is equivalent to less than 18 % of the pore volume.
One should keep in mind that this result was achieved with a single vertical injection
well and a reservoir height of 100 m. More wells, different setups, or different reser-
voir geometry obviously can lead to a more efficient or less efficient utilisation of the
available pore space.
• The concept of calculating storage capacity for a model domain with a horizontal cap
and one injection well is a conservative approach since structural-, solubility-, and
chemical trapping mechanisms, which are not considered here, would lead to further
increase in storage capacity. However, these mechanisms mainly contribute to total
trapping at later times.
• High storage capacity of a given reservoir is typically achieved for low injectivities.
This is in conflict with an economically feasible performance, and compromises have
to be made in engineering practise.
6 Sensitivity Analysis
In the previous chapters, a sound understanding of the forces and processes influencing
the plume evolution and the resulting storage capacity has been gained. Detailed storage-
capacity estimates have been calculated by numerical simulations for reservoirs employing
typical reservoir parameter setups. The question now arises, as to the influence of individ-
ual parameters like porosity, permeability, residual saturations, etc. on the model results,
considering their full parameter range. The scope of this chapter is to determine the influ-
ence of various reservoir parameters on the model results with respect to storage capacity
and risk and to order qualitatively the parameters depending on their influence. This is a
prerequisite for the following Chapter 7, where a risk analysis is performed. Since in a risk
analysis the computational cost is very high, only parameters with the strongest influence
can be chosen to be independent, whereas the others (of less influence) are calculated by
functional dependencies.
Sensitivity Analysis methods can be grouped into three classes (Saltelli et al., 2000): “screen-
ing”, “local” and “global”.
87
88 6.1 Discussion of Sensitivity Analysis
measures, i.e. they provide a ranking of input factors, but do not provide quantitative
sensitivity measures.
• Local sensitivity analysis methods give local quantitative measures of input factor
influence. In other words, they often compute partial derivatives at one point in
multidimensional parameter space. Thus, the sensitivity measure is only valid locally,
and especially for highly non-linear models, this produces a limited applicability. With
respect to the computational cost of the evaluation, they range between screening and
global sensitivity analysis methods.
• Global sensitivity analysis methods apportion the output uncertainty to the uncer-
tainty in the input factors. This approach usually varies all input factors simulta-
neously and sensitivity measures are given over the entire range of each input factor
distribution. The computational cost for a global sensitivity analysis is high; more
than a hundred model evaluations are necessary per input factors.
For all the sensitivity analysis methods considered, there is a clear trade-off between com-
putational cost and possible information to be gained. With respect to the numerical ex-
periments undertaken in this section for an analysis of CO2 injection processes, the com-
putational cost can be considered very high. The complexity of the model together with
the highly non-linear fluid properties (cf. Section 2.2.2) and constitutive relations (cf. Sec-
tions A.1.3 and A.2.4) impose a highly non-linear model result behaviour. On the one hand,
this non-linear behaviour would require choosing a global sensitivity analysis method; on
the other hand, it precludes it (due to the high computational cost). A hint for solving
the dilemma is given by Campolongo et al. (1999). They propose performing a screening
analysis based on the so-called Morris Method. The Morris Method can be considered a
“global screening analysis”. It can be regarded “screening” because one input factor out of
many is modified in each model evaluation. This is known as a one-at-a-time design. It is
“global” because the methodology randomly selects points in the input factor space. Thus,
a qualitative input factor ranking is produced with respect to the entire input factor space.
One severe drawback, however, remains: no quantitative sensitivity measures are produced.
xi have been normalised to a range of variation of unity (xi ∪ (0, 1)). Thus, the p-level values
1 2
of xi can be written as {0, p−1 , p−1 , . . . , 1}. To obtain the input factor for an experiment x,
one individual input factor xi is varied by ∆, a predetermined level value, i.e. a multiple of
1
p−1
. The outcome of a (numerical) experiment is written as y = y(x1 , . . . , xk ).
The sensitivity measure proposed by Morris is based on so-called “Elementary Effects” (EE).
For the ith input factor EEi is calculated to
y(x1 , . . . , xi-1 , xi + ∆, xi+1 , . . . , xk ) − y(x1 , . . . , xk )
EEi (x1 , . . . , xk ) = . (6.1)
∆
Note that EEi is a local measure of sensitivity at point x. To obtain a global measure of
sensitivity, several randomly selected points xr in the input parameter space are evaluated,
where r is the number of local measures. In his original study, Morris proposes two measures
of sensitivity for each input factor:
A high value of the mean of the distribution of the Elementary Effects (µ(xi )) (compared
to µ of other input factors) indicates an overall high influence on the outcome of the ex-
periment by input factor i. A high value of the standard deviation of the distribution of
the Elementary Effects (σ(xi )) (compared to σ of other input factors) indicates interaction
effects of the considered factor with other factors, or a non-linear effect on the outcome of the
experiment. To put it precisely, (σ(xi )) does not provide estimates of individual interaction
among different input factors, but it estimates whether any significant interaction exists.
Campolongo and Braddock (1999) realised that, when Elementary Effects EEi have opposite
signs, they may cancel each other out. They proposed a modified measure for the mean of
the distribution of the Elementary Effects by only considering absolute values:
To save computational costs, so-called “trajectories” are constructed. This is the selection
of the point x as a starting point to estimate the Elementary Effects of an input factor as
the endpoint of the preceding estimation:
y(x1 + ∆, x2 ) − y(x1 , x2 )
EE11 = (6.5)
∆
y(x1 + ∆, x2 + ∆) − y(x1 + ∆, x2 )
EE21 = (6.6)
∆
Thus, xr is only selected randomly for each trajectory, not for every experiment. The overall
cost of the method can be calculated as r · (k+1) where r ∼ 4–10.
The Morris Method in the version explained above has been tested by Campolongo et al.
(2004); the authors state that the method is efficient in identifying irrelevant factors (i.e.
low µ∗ and σ). Cropp and Braddock (2002) demonstrated that analyses even using quite a
90 6.2 Parameters
small number of trajectories (r) at coarse resolutions (low number of levels p) can provide
very good estimates of model sensitivity to parameter perturbations and interactions of up
to three parameters. They further state that the method is a powerful and efficient tool for
global sensitivity analysis.
The software package SimLab (Joint Research Centre of the European Commission, 2004)
is a development framework for uncertainty and sensitivity analysis and is used to generate
the input parameter sets and is used for the analysis of the model output results.
6.2 Parameters
The reservoir parameters investigated in the sensitivity analysis are given in Table 6.1 along
with the minimum and maximum values allowed and the literature source. For parameters
kh , φ, dip angle, S, dT/dz, and D, the 5th (Min) and 95th (Max) percentile values of the
NPC-database statistics (cf. Table 3.1) are used to limit the range. The maximum capillary
entry pressure value pd and the maximum residual water-rich phase saturation Sw,r are the
maximum values measured by Bennion and Bachu (2008). The minimum CO2 -injection
temperature (TCO2 ) is selected as the critical temperature. The injection interval ranges
from injection only in the lowest 10% of the reservoir height to a CO2 injection over the
entire reservoir height (100%). Other values are best-guess values.
Table 6.1: Parameters investigated in the sensitivity analysis, parameter range, and literature
source.
When these 15 input parameters are considered in a Morris Method-type sensitivity inves-
tigation, the number of necessary model runs using four trajectories is 64. The parameter
distribution within the range is assumed to be uniform. To avoid an influence of different
6 Sensitivity Analysis 91
input parameter ranges on the sensitivity measures, all parameter ranges are normalised to
a range of zero to one. Consequently, the results of the numerical simulations have also been
normalised to a range of zero to one to calculate µ∗ and σ.
6.3.2 Results
In Figure 6.2 (left), the centres of gravity of the CO2 plumes of the individual model runs
are shown as diamonds. A large variability of resulting CO2 mass distributions within
the reservoir is observed. There are combinations of input parameters, leading to a CO2
plume spreading far into the reservoir (large centre of gravity x value), others lead to a CO2
plume closely distributed around the injection boundary (small centre of gravity x value).
Consequently, the information value of the derived sensitivity measures are representative
92 6.3 Numerical Investigations
Length
Injection
Interval
Dip angle α
for a broad range of flow and transport conditions in the reservoir. The simulations are
stopped before the CO2 plume reaches the eastern boundary. This is to prevent unwanted
influence of boundary conditions on the results.
101
gravitational
forces
increase
100
100
viscous
forces
80 increase
Centre of Gravity y[%]
-1
10
60
Gr [-]
10-2
Injection Interval
40
-3
10
20
viscous capillary
10-4 forces forces
0 increase increase
10-5 -5
0 2 4 6 8 10 10 10-4 10-3 10-2 10-1 100 101
Centre of Gravity x[%]
Ca [-]
Figure 6.2: Centre of gravity of the CO2 plume given in percent of the length of the domain
(2000 m) versus percent of the height of the domain (H) (left). Capillary Number Ca
versus Gravity Number Gr (right). Both figures are shown after a total of 22 tCO2 /m
is injected. Each diamond indicates an individual model run.
Figure 6.2 (right) shows Capillary Number Ca versus Gravity Number Gr. The dimensionless
numbers have been calculated by using CO2 -mass weighted averages over all discretisation
cells for CO2 density, viscosity, capillary pressure and velocity. Furthermore, k is repre-
sented by horizontal permeability. Characteristic length lcr is represented by the height of
the domain D. The actual front widths have not been considered, which might lead to an
error in the range of less than one order of magnitude in Ca. The range of variation for Ca
and Gr is three orders of magnitude, mainly due to a variation in absolute permeability kh ,
changing by three orders of magnitude. A vague trend can be identified from high Gr / high
6 Sensitivity Analysis 93
Ca to low Gr / low Ca. This behaviour is already observed in the analytical experiments
in the dimensional analysis for varying characteristic velocities vcr (cf. Figures 4.2 and 4.9
(right)). This relation of dimensionless numbers (balance of forces) represents a snapshot
after 22 tCO2 /m have been injected. The numerical experiments here are rather dominated
by viscous forces. For model runs with conditions leading to dimensionless numbers close
to the equilibrium of forces (Gr and Ca close to one), gravitational and capillary forces gain
importance. Hence, parameter sensitivities can only be derived and interpreted for this state
of the system; i.e. parameter sensitivities might be different for a reservoir that is dominated
by gravitational or capillary forces.
To perform a sensitivity analysis, the questions of interest need to be clearly defined. The
focus here is to define questions which are important to answer in storage capacity studies and
in risk assessment studies. Hence, the parameter sensitivities with respect to the answers to
these questions enables the identification of parameter importance in these types of studies.
The following questions are addressed:
• Which fraction of the total injected CO2 is in free-phase (not dissolved in brine) in the
topmost region of the reservoir at a given point in time?
The fraction of free-phase CO2 mass in the topmost region of the reservoir is a measure
of the strength of the buoyancy forces. As shown in Section 5, strong buoyancy forces
reduce storage capacity. The mass fraction is also a measure of risk, since only the
free-phase CO2 in the topmost regions could potentially leak if a leakage pathway is
present.
• What is the mass ratio between CO2 occurring in free phase and CO2 dissolved in
brine at a given point in time?
Here, the same reasoning holds with respect to risk as in the previous bullet, only the
free-phase CO2 could potentially leak if a leakage pathway is present. Carbon dioxide
dissolved in brine can be considered as safely stored.
To answer the first question, an observation point needs to be defined. A point at 20 metres
distance from the western boundary and 20 metres above the bottom of the domain is se-
lected. It is necessary here to select a point rather close to the western injection boundary
94 6.3 Numerical Investigations
due to the broad parameter range considered. On the one hand, there are cases with pa-
rameter combinations, which lead to strong gravity segregation, i.e. quick upward movement
of the CO2 plume and fast lateral spreading below the caprock. On the other hand, there
are cases with parameter combinations having e.g. small injection intervals and weak grav-
ity segregation. In these cases, the CO2 plume does not reach the caprock at all and CO2
spreads only a small distance into the reservoir. Hence, the observation point selected is the
most distant point from the injection boundary where the CO2 arrival is observed for all
cases.
To answer the second and third questions, a time needs to be defined. Here, the sensitivity
measures are compared at the point in time when the same amount of CO2 is stored in the
reservoirs, i.e. 22 tCO2 /m.
The input parameter sensitivity measures µ∗ and σ (cf. Section 6.1.1) for the questions stated
are given in Figures 6.3 and 6.4.
1 K 1Permeability A Anisotropy
A Anisotropy B Depth
P Porosity D Dip angle
0.8 0.8Height
H E Pd
D Dip angle H Height
S Salinity I InjTempCO2
0.6 0.6 W
E Pd K Permeability
σ
P B MR
Depth T Temperature Gradient
SI
S
L
D
0E
AINH
R W D
KP
0Injection
L Interval W Injection Interval
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
µ µ
* *
Figure 6.3: Input parameter sensitivity measures µ∗ (e.g. Equation 6.4) and σ (e.g. Equation 6.3)
with respect to CO2 arrival time at an observation point (left) and with respect to
the total CO2 mass (free phase) in the upper 20 % of the reservoir (right). Low µ∗
indicates little overall influence on model answers. Low σ indicates parameters with
little interaction of the considered factor with others.
Figure 6.3 (left) shows sensitivity measures with respect to the arrival time of the CO2 plume
at the observation point. The arrival time at this point is sensitive to permeability, injection
interval and injection rate, geothermal gradient and reservoir depth. Sensitivity to the first
two parameters can be attributed to the developing overpressure, as will be shown later.
The sensitivity of the injection interval is due to the close position of the observation point
to the injection boundary. Geothermal gradient and reservoir depth are sensitive parameters
since they account for the density of the CO2 , and therefore for the buoyancy forces. Other
input parameters are negligible. All parameters mentioned show a comparable overall effect
6 Sensitivity Analysis 95
Figure 6.3 (right) shows sensitivity measures with respect to the fraction of the total injected
CO2 in free phase located in the topmost 20 % pore space of the reservoir after 22 tCO2 /m
have been injected. The injection interval has the highest overall effect on the model answers.
This can be explained by a simple example. If the injection interval is 1.0, at least 20 % of
the CO2 mass are located in the topmost 20 % of the pore space of the reservoir just after the
injection start. If the injection interval is less then 0.8, no CO2 is located in the topmost 20 %
of the pore space of the reservoir just after the injection start. The amount of CO2 located in
the topmost region after some time thus depends not only on the effects of buoyancy forces;
especially at early times, it can largely be influenced by the injection interval. The strength
of the buoyancy forces is largely influenced by the density difference between the water- and
CO2 -rich phases, these depend on the depth, the geothermal gradient, and the permeability
anisotropy. These parameters are shown to have a large effect on model answers. The model
answers are also sensitive to the height of the reservoir, since it accounts for the distance the
CO2 has to migrate in order to reach the topmost 20 % of the pore space of the reservoir.
Some influence is also found for the capillary entry pressure and the residual CO2 -rich phase
saturation. The residual CO2 -rich phase saturation parameter accounts for the amount of
CO2 which is residually trapped, and therefore reduces the amount of CO2 which may reach
the topmost 20 % of the pore space of the reservoir.
1 K 1Permeability A Anisotropy
A Anisotropy B Depth
P Porosity K D Dip angle
0.8 0.8Height
H E Pd
D Dip angle H Height
S Salinity I InjTempCO2
0.6 0.6
E Pd K Permeability
σ
Figure 6.4: Input parameter sensitivity with respect to the mass fraction of CO2 occurring in free
phase in relation to the CO2 dissolved in brine (left) and with respect to maximum
overpressure developing below the caprock (right). Low µ∗ indicates little overall
influence on model answers. Low σ indicates parameters with little interaction of the
considered factor with others.
Figure 6.4 (left) shows sensitivity measures with respect to the mass ratio of CO2 occurring
in free phase and CO2 dissolved in brine when 22 tCO2 /m are injected. The mass ratio de-
96 6.3 Numerical Investigations
scribed is not dominated by a single or few parameters. Rather five to ten parameters have
some effect. Among these influencing parameters are permeability, depth, residual CO2 -rich
phase saturation, geothermal gradient, salinity, sorting factor, residual water-rich phase sat-
uration injection interval, and height of the reservoir. It is surprising that the mass injection
rate has very little effect.
Figure 6.4 (right) shows sensitivity measures with respect to the maximum overpressure
developing in the reservoir due to CO2 injection. The maximum overpressure developing
is sensitive to the absolute permeability, the CO2 injection rate and the height of the do-
main, in this order of importance. Other input parameter influences are negligible here.
The absolute permeability shows the highest overall influence (high µ∗ ) on the developing
overpressure and has either a strong non-linear influence on the overpressure or interaction
effects with other parameters (high σ), possibly CO2 injection rate. The developing over-
pressure is sensitive to the height of the domain because it is comparable to the available
pore space the CO2 is injected into. One would expect some sensitivity of overpressure with
respect to the CO2 injection temperature or depth, since these parameters influence CO2
density and thus CO2 compressibility, but they seem to be negligible when compared to the
influence of permeability, injection rate and domain height. This information is of practical
use when the overpressure below a caprock is to be limited. Then it is most advisable to
keep injection rates low and search for highly permeable reservoirs, which is what one would
expect. Varying the CO2 injection temperature and the injection interval does not have an
effect.
Assuming that these questions are representative and allow for conclusions on parameter
importance in storage capacity estimation and in risk assessment studies, a qualitative pa-
rameter ranking can be derived, as shown in Table 6.2. This is done by averaging the overall
parameter effect in the four individual parameter sensitivity rankings.
According to Table 6.2, parameters like horizontal permeability, injection interval, reservoir
depth, reservoir height, and geothermal gradient have the largest influence on the model
results with respect to the questions stated here. On contrast, parameters like the salinity,
residual water-rich phase saturation, entry pressure, porosity, sorting factor, CO2 injection
temperature and the dip angle can be neglected in further investigations in the fields of
storage capacity estimation and risk assessment. However, this is only valid for the chosen
setup. It was shown that the setup chosen here is strongly influenced by viscous forces. In
other setups, which might be strongly influenced by gravitational or capillary forces, the
parameters identified here as having little influence might become important. This has to
be discussed when neglecting parameters.
6 Sensitivity Analysis 97
increasing importance
7 SCO2 ,r Residual CO2 -rich phase saturation
8 AnIso Permeability anisotropy
9 S Salinity
10 Sw,r Residual water-rich phase saturation
11 pd Entry pressure
12 φ Porosity
13 λBC Sorting factor
14 TCO2 CO2 injection temperature
15 α Dip angle
Table 6.2: Qualitative ranking of parameter effect as the average of individual overall parameter
effects for the four stated questions of interest.
98 6.4 Summary and Conclusion
• Parameters like the salinity, residual water-rich phase saturation, entry pressure, poros-
ity, sorting factor, CO2 injection temperature and the dip angle have little influence
on the model results with respect to the questions stated here and with respect to the
setup chosen.
• For other setups than the one considered here, parameters identified here as having
little influence might become important. For example in two-dimensional space, the
increase of pore space with increasing distance r from the injection boundary is lin-
ear (H · r). In three-dimensional space, an exponential increase of pore space with
increasing distance from the injection boundary is observed (H · π · r2 ). In the latter
case, this results in a rapid decrease of total velocity and therefore a rapid decrease of
viscous forces which dominate here. With an increasing influence of buoyancy and/or
capillary forces, other parameters, which have been identified as being negligible here,
could gain importance.
7 Risk Analysis∗
In the previous Chapter 6, a sensitivity analysis was performed to order reservoir parameters
qualitatively, depending on their influence on storage capacity and risk. In the following risk
analysis, only the parameters with the strongest influence are chosen to be independent,
whereas the others (with less influence) are calculated by functional dependencies. The rea-
son for this simplification is the high computational cost that is generated by considering
all parameters to be independent. The risk analysis investigates a scenario where CO2 is
leaking out of a storage reservoir through a pre-existing leaky well. Risk is expressed in this
case as the likelihood of this leakage to occur times the leaked CO2 mass at a given point
in time and for a given distance between the leaky well and the injection well. Numerous
simulation experiments are run, employing individual reservoir parameter setups. To set
up an individual simulation experiment, three out of four independent parameters are ran-
domly sampled from the database parameter distributions given in Chapter 3. The fourth
independent parameter is randomly sampled from a parameter distribution calculated by a
theoretical model. Thus, the risk analysis is performed with respect to the database reser-
voir parameter distributions, by considering a theoretical model, and by selected literature
sources on functional dependencies of the dependent parameters.
While the concept of risk developed and applied in this study is straightforward, it is intended
to provide
• a systematic framework for engineers for determining the ideal properties for a subsur-
face CO2 storage reservoir,
• assistance in rating storage sites in terms of risk depending on their properties. This
knowledge can assist experts when utilising more comprehensive screening and ranking
frameworks,
• estimates of risk of a site which can be compared to the risk of other sites, possibly
leading to the decision where to conduct further investigations,
• a selection of the most appropriate location for an injection well based on risk estima-
tions depending on the number and distribution of surrounding leaky wells.
∗
This Chapter is submitted to Energy Conversion and Management with the title “Kopp, A., Binning,
P.J., Johannsen, K., Class, H. and Helmig, R., Risk Analysis for Leakage through Abandoned Wells in
Geological CO2 Storage” by December, 17th 2008.
99
100 7.1 Discussion of Risk
two time scales may be of interest, and these are directly related to the environmental im-
pacts that might occur from leakage: local environmental effects or global effects. The scope
of the risk assessment and the associated time scales are defined dependent on these impacts
(IPCC, 2005). The medium to long timescale (∼20 years or longer) considers global effects,
that is the return of the CO2 back into the atmosphere. The short time scale considers local
environmental impacts, leading to risk directly associated with exposure to leaked CO2 ; that
is local health, safety and environmental hazards. For example, Pruess (2008b) suggests a
potential for self-enhancement of leakage rates, leading to a so-called “pneumatic eruption”;
although the author says it is unlikely. Pacala (2003) states that local risks (together with
economic considerations) are likely to constrain allowable leakage rates more tightly than
global impacts. Hence, the focus of the risk assessment method developed here is on the
short time scale, i.e. up to ∼20 years.
Relating to the first question, there are many ways in which a CO2 storage system can
fail. For example the system can fail in providing the necessary injectivity with respect
to the CO2 mass delivered by the surface installations (power plant, pipe network, process
engineering devices etc.). There are a number of engineering failures associated with the
design of the installations and the injection process, for example well corrosion, untight well
plugs, formation clogging due to halite precipitation in the close vicinity of the injection
102 7.1 Discussion of Risk
well. There might be failure associated with the design of the CO2 monitoring devices to
measure CO2 plume evolution and possible CO2 leakage. The system can also fail to provide
the necessary capacity to store the intended CO2 production. Finally, the system can fail
after injection due to CO2 leakage, e.g. through the injection well.
The risk study conducted here focuses on CO2 leakage out of the storage reservoir. IPCC
(2005) summarised number of possible leakage pathways: (i) through the pore system in low-
permeability caprocks such as shales, if the capillary entry pressure at which CO2 may enter
the caprock is exceeded; (ii) through openings in the caprock or fractures and faults; (iii)
through anthropomorphic pathways, such as poorly completed or abandoned pre-existing
wells. The focus is here on leakage through poorly completed or abandoned pre-existing
wells. Abandoned wells have been identified as one of the most probable leakage pathways
for CO2 (Benson (2005), Gasda et al. (2004)). A pre-existing (leaky) well is assumed to exist
at some distance from the injection well. It has a constant diameter of 1 m and is screened
over the full reservoir thickness. Other (presumably smaller) well diameters can easily be
considered, as will be shown later. The leakage pathway in the abandoned well is not of
importance, as will be explained in the definition of the consequences of failure.
The likelihood of failure can be assessed by considering the potential range of reservoir
parameters like porosity, geothermal gradient, depth, and anisotropy. The statistical char-
acteristics of reservoir parameters were determined from a database including information
from over 1200 reservoirs or by a model (permeability anisotropy). The parameter space of
reservoir properties is randomly sampled and simulations are conducted to assess the dis-
tribution of CO2 in the subsurface. For each simulation, failure is assessed by examining
whether CO2 has spread to a given radius within a given time. If CO2 has spread beyond
the leaky well radius, this case has failed. By simulation of many such cases, a likelihood of
failure can be given by relating the number of cases that failed nf (dependent on distance r
from the injection well to the leaky well and time t) to the total number of cases simulated
(N).
To determine the consequences of failure, the damage is defined as the amount of CO2 mass
that has spread beyond the leaky well distance (compare to Figure 7.1). The assumption
is that damage (Di ) measures the mass that can potentially leak through the pre-existing
well out of the system. This means that in this study the leakage process itself is not
modelled. Consequently it is not necessary to be more specific about the leakage pathway
in the abandoned well, as given in Gasda et al. (2004).
IPCC (2005) outline a “framework for assessing environmental risks”. They define two
categories of environmental impacts that might occur from geological CO2 storage: local
environmental effects and global effects. Global effects arise from CO2 leakage to the atmo-
sphere and relate to the uncertainty in the effectivity of CO2 storage. Local effects arise from
the direct exposure of surrounding plant and life species to CO2 or other dangers associated
with the leakage in the close surrounding of the site. This study partially addresses both
impact categories. The CO2 mass that has spread beyond the leaky well can potentially leak
through the well. However, the impact of this CO2 is not explicitly considered. The leaked
7 Risk Analysis 103
CO2-front Leakage
mCO2
rleaky well
r
Figure 7.1: Sketch of the radially symmetric model domain and definition of the leakage that
occurred after time t at a leaky well in distance rleaky well from the injection well. CO2
is injected into the centre boundary and spreads laterally in the reservoir. Buoyancy
forces drive the CO2 upwards where it accumulates below the caprock. Consequently
the CO2 plume spreads faster at the top than at the bottom.
CO2 could possibly leak (i) to the next shallower geological unit where it is safely stored, or
(ii) leak to a fresh-water bearing aquifer and pollute drinking water, or (iii) leak back to the
atmosphere. In cases (ii) and (iii) it is important to determine the rate at which the CO2
leaks. If it leaks at considerable rate, local health, safety and environmental hazards may
arise, yielding risk that is much higher and may not be proportional to the amount of CO2
that has leaked (IPCC, 2005).
Since in this study a simple reservoir geometry is used together with a random sampling of
statistical parameters distributions, this risk assessment does not refer to any specific site.
Risk, as it is defined here, refers rather to a statistical leakage probability that might occur
at a pre-existing well. Summarising, risk is defined by multiplying the likelihood of failure
with the consequence of failure, which is expressed in Equation 7.1.
N
nf X Di
Risk [kg] = · , (7.1)
N
|{z} i=1
nf
likelihood of failure | {z }
consequence of failure
where nf is the number of cases that failed [-] (dependent on distance r [m] from the injection
well to the leaky well and time t [d]), N is the total number of cases simulated [-], and Di is
the damage [kg] of case i at time t.
Note that since a radially symmetric domain is assumed with a leaky well diameter of 1 m
(as indicated in Figure 7.1), the damage Di is calculated by simply multiplying it with
1m
2πr
, which represents the fraction of a full circle affected by the leaky well, where r is the
distance of the leaky well from the injection point. To consider smaller leaky well diameters,
104 7.2 Parameters
the resulting damage and risk values presented in the following just need to be multiplied by
the leaky well radius (since results are given for a unit leaky well diameter, i.e. 1 m). Risk
has the unit of mass (kg). Since leakage is not explicitly modeled (no leaky well is included
in the simulations), a single model run can be used to calculate damage and risk for many
leaky well radii r.
The concept outlined above, calculates risk as an average over the entire parameter space.
This becomes apparent because nf cancel out of Equation 7.1 and what is left is average risk
PN Di
i=1 N . The average risk, using the statistical distribution within the parameter space, can
be used to draw some important conclusions (as given in Sections 7.3.3 and 7.5), but does
not determine the risk of an individual case.
A procedure to calculate risk for an individual case, could be the following: Assume that each
simulation has four independent input parameters, p1 . . . p4 , among the parameter distribu-
tions. When fixing p2 , p3 , and p4 , the probability of failure P1 due to variation of parameter
p1 can be determined by varying p1 and finding p∗1 that partitions between failure and no
failure. The probability P1 is then determined from the probability distribution to be the
probability that p1 > p∗1 , where p1 > p∗1 defines the range of parameters p1 leading to failure.
The total failure probability P of this case can then be calculated by P = P1 · P2 · P3 · P4 .
The risk associated with this case would be the damage produced by this case times the
failure probability, i.e. riski = Di · P1 · P2 · P3 · P4 . This individual risk is of course dependent
on the distance r from the injection well to the leaky well and the time t. This study does
not consider the risk of an individual case any further, but instead focuses on the average
risk and the conclusions that can be drawn from it.
7.2 Parameters
To define parameter input sets for simulations, four independent parameter distributions are
randomly sampled. These four parameters are porosity φ, depth D of the reservoir below
surface, average geothermal gradient dT/dz at the site, and the anisotropy between vertical
and horizontal absolute permeability AnIso. These are called primary parameters. All other
parameters, called secondary parameters, are functions of primary parameters or constants.
The primary parameters have been selected based on previous sensitivity investigations
(cf. Section 6.3.2) where they have shown to be the parameters which are most influential
on CO2 plume evolution behaviour in a reservoir† . When comparing the qualitative ranking
of parameter effects (cf. Table 6.2) determined in the sensitivity analysis, primary param-
eters rank on third, fifth, eighth and twelfth position. Other parameters are not selected
as primary parameters for three reasons, (i) there is a lack of data, (ii) there is a mutual
interrelation with primary parameters, (iii) computational effort would be too high. Detailed
†
Kopp, A., Class, H. and Helmig, R., Sensitivity Analysis of CO2 Injection Processes in Brine Aquifers,
presentation given at European Geosciences Union (EGU) General Assembly, April 18th , Vienna, Austria,
2007
7 Risk Analysis 105
∗
%b is constant for the calculation of initial and boundary conditions (1160 [kg/m3 ]).
Table 7.1: Definition of primary and secondary model input parameters, dependencies, and
sources.
7 Risk Analysis 107
To calculate statistical characteristics of the 4th primary parameter, the absolute perme-
ability anisotropy (AnIso), a layered reservoir is considered. Layers have varying thickness
and an alternating constant absolute permeability of 3.24 · 10−15 m2 and 1451.85 · 10−15 m2 .
These permeabilities reflect the 5th and 95th percentiles of the absolute permeability distri-
bution in the NPC database. When estimating upscaled or effective permeability for flow
in a layered reservoir, the direction of the flow is of importance. The harmonic mean HM
of the permeabilities is used to estimate effective permeability for flow perpendicular to the
layers:
P
mi
HM = Pi mi , (7.2)
i ki
where i is the layer index, mi [m] is the thickness of the layer having permeability ki .
The arithmetic mean AM of the permeabilities is used to estimate effective permeability for
flow parallel to the layers,
P
i mi · ki
AM = P . (7.3)
i mi
The anisotropy of effective flow permeability is calculated as the harmonic mean divided
by the arithmetic mean. When varying the number and the thickness of the layers, the
distribution as given in Figure 7.2 is the result.
Statistical characteristics of the distribution are: minimum = 0.0089, maximum = 1.0,
arithmetic mean = 0.0283, median = 0.0118, 5th percentile = 0.00891 and 95th percentile =
0.0842.
Section 7.2.3 gives the procedure of sampling the primary parameter distributions in order
to define a simulation case.
relative
frequency
0,40
0,30
0,20
0,10
0,00
-2,0 -1,5 -1,0 -0,5
Ln (Anisotropy) [-]
Figure 7.2: Histogram data show relative frequency of the natural logarithm of anisotropy between
vertical and horizontal intrinsic permeability [-] derived from a anisotropy model.
Line indicates a normal distribution having the same statistical characteristics as the
histogram data set.
Correlation coefficients of functions and database values are rather poor. However, the
Kozeny-Carman function given in Equation 7.4 with an average grain diameter of 35 µm fits
the interpolated NPC database values reasonably well. The absolute horizontal permeability
is therefore calculated from porosity by using:
1 φ3
kh = D50 , (7.4)
150 (1 − φ)2
where D50 is average grain diameter [µm]. Hence, absolute permeability, which is the pa-
rameter with the highest effect in the sensitivity analysis (cf. Table 6.2), is considered by
a dependency on porosity. Porosity was preferred as a primary parameter, since porosity
measurements are more reliable than permeability measurements.
Injection Interval (II), CO2 injection rate qCO2 , and reservoir height are assumed
to be constant in this analysis, although they showed a high effect in the sensitivity anal-
ysis (rank two, four and six in Table 6.2). The reason for not choosing them as primary
parameters is, that the injection interval and the CO2 injection rate can easily be modified
when drilling and screening one or more injection wells. They are not a reservoir property.
The reservoir height is a reservoir property, here a lack of data prevents consideration as a
primary parameter.
Residual water- and CO2 -rich phase saturations (Sw,r and SCO2,r ) are defined as con-
stant since no suitable correlation function is known to the author. Holtz (2002) gives a
7 Risk Analysis 109
1,E-09
1,E-10
1,E-11
1,E-12
1,E-13
1,E-14
1,E-15
1,E-16
5 10 15 20 25 30 35 40 45
Porosity [%]
Figure 7.3: Correlation functions between porosity and absolute permeability and NPC database
reservoir values. Given are correlations by a Kozeny-Carman model (Scheidegger
(1960)) for different average grain diameters, correlations by Pape et al. (1999) for
average sandstones and Rotliegend sandstone, a correlation given by Holtz (2002),
and an interpolation of the NPC database values (k [m2 ]=2.046 · 10−11 φ3.6555 ).
correlation function for the residual water-rich phase saturation dependent on absolute per-
meability and porosity for a specific site. However, permeability and porosity ranges and
distributions leading to the functionality given, are not identical to the data considered here.
This results in unrealistic values for residual water-rich phase saturation Sw,r .
Relative permeability kr is defined by a Brooks & Corey model (Brooks and Corey (1964)).
Since Brooks & Corey input parameters Sw,r , SCO2,r together with the sorting factor λBC are
constant, the relation is identical in all simulations. This is a simplification which is neces-
sary due to lack of better knowledge (data). Only very few measured relative permeability
relations of CO2 -brine systems are known, e.g. Bennion and Bachu (2005). But these sparse
data do not allow the definition of dependencies on e.g. primary parameters. Relative perme-
ability has a high influence on plume evolution behaviour, as will be discussed in Section 7.4.
water-rich and the CO2 -rich phase, water-rock contact angle, and water-rich phase satura-
tion. A measured capillary pressure-saturation relation and the Leverett J-function (Lake,
1989) are used to determine a distribution of capillary pressure relations as a function of the
primary parameters. The measured data are from Plug and Bruining (2007) experiment no.
12 where the sand column has a permeability of approx. 2 · 10−10 m2 and a porosity of 0.37.
The experiment was conducted at a pressure of 8.5 MPa and at 300.15 K. The author also
gives values for interfacial tension σ and contact angle Θ. These data allow to use a Leverett
J-function to normalise the capillary pressure relation. This relation is representative for the
investigated rock-type. The dimensionless Leverett J-function is given in Equation 7.5:
s
pc (Sw ) kh
J(Sw ) = . (7.5)
σ · cosΘ φ
The assumption, underlying Equation 7.5 is, that theqporous medium can be modelled as a
bundle of non-connecting capillaries. In this model, kφh [m] is proportional to the typical
pore throat radius. The Leverett J-function is used to scale the given capillary pressure rela-
tion to other permeability, porosity, temperature and pressure (since the interfacial tension
σ changes with pressure and temperature). This is done by equating the J-function for two
cases and solving for the capillary pressure of interest (cf. Equation 7.6).
s
kh,1 φ2 σ2 cosΘ2
pc,2 (Sw ) = pc,1 (Sw ) . (7.6)
kh,2 φ1 σ1 cosΘ1
Subscript “1” indicates reference capillary pressure relation and rock properties (here ex-
periment 12 in (Plug and Bruining, 2007)), subscript “2” the capillary pressure relation for
rock properties of interest. An equation for the interfacial tension σ for water and CO2
systems, as a function of pressure and temperature, is given by (Kvamme et al., 2007) (see
Figure 2.10). Examples of resulting capillary pressure relations are shown in Figure 7.4.
0.4
0.2
0.1
0
0.1
0 0.2 y [-]
0
0.2 0.3 ro s it
0.4
Water S 0.6 0.4 Po
aturation 0.8 1
[-]
Figure 7.4: Capillary pressure dependence on water-rich phase saturation and porosity. Here
σ1 = σ2 , i.e. relations are at constant temperature and pressure.
and results will be discussed with respect to this simplification in Section 7.4.
CO2 fluid properties (density and viscosity) are calculated dependent on temperature and
pressure. Since in the risk analysis reservoir depth D is a primary paramater and its range
includes depths where CO2 occurs in gaseous, liquid and supercritical state of aggregation,
the calculation of CO2 fluid properties needs to cover a broad range of possible pressures
and temperatures.
Other secondary parameters, like brine fluid properties, solubility of CO2 in brine etc.
are calculated by functions given in literature (as cited in Table 7.1) or are defined as con-
stant values since results are insensitive to variation of these parameters.
All parameters given in the NPC database are tested for mutual interrelations (cf. Sec-
tion 3.3). Apart from the considered dependencies (between porosity and absolute perme-
ability and between depth and temperature, which is the geothermal gradient) no significant
correlations have been found.
is obtained. Primary parameters are constant throughout a simulation. Following the gen-
eration of four random primary parameters, the secondary parameters are calculated. All
secondary parameters, which are dependent on the unknowns in the simulation, i.e. pressure,
temperature, saturation or CO2 mass fraction dissolved in brine, are continuously updated
during model simulations. This means, for example, that the capillary pressure relation could
be different in each box (cf. Section A.4) and that it changes with time. Other secondary
parameters are fixed during the simulation.
7.3.3 Results
In Figure 7.5 typical results are shown for two cases. The amount of injected CO2 is in
both cases 20 Mt CO2 . Randomly sampled primary parameters for Case 11 are φ11 = 0.19,
dT/dz11 = 0.020 ◦ C/m, D11 = 1920.2 m, AnIso11 = 0.009 and for Case 48 are φ48 = 0.26,
dT/dz48 = 0.053 ◦ C/m, D48 = 640.1 m, AnIso11 = 0.023. These settings lead to CO2 density
around 790 kg/m3 in Case 11 and 330 kg/m3 in Case 48. In Case 11 no leakage occurs up to
the time represented by the figure, whereas for Case 48 leakage occurs.
Figure 7.5: Slice of the radially symmetric domain showing CO2 -rich phase saturation after
20 years of continuous injection for Case 11 (left) and Case 48 (right). Shown is
only the centre (inner) part of the domain where reservoir height is exaggerated. No
leakage can occur in Case 11, whereas for Case 48 CO2 has spread beyond leaky well
distance and potential leakage is high.
Given the definition of risk in Equation 7.1, the two terms of the equation are examined
separately before looking at the product of both terms, yielding risk.
In Figure 7.6 the likelihood of failure is shown. Each point on the surface resembles the
combined results of all simulations conducted. Range of the likelihood of failure is between
zero (no case has failed) and one (all cases have failed). Initially no CO2 is injected and
consequently no failure has occurred. After 7000 days all cases below rleaky well = 1100 m
have failed and the likelihood of failure is equal to one. For distances larger than rleaky well =
1100 m the likelihood of failure is smaller, since less cases fail with increasing distance. For
rleaky well = 2000 m only 20 % of all cases fail after 7000 days of CO2 injection. Generally the
likelihood of failure increases with time and decreases with distance rleaky well .
In Figure 7.7, the consequence of failure is shown. Each point on the surface resembles the
combined results of all simulations conducted. Range of the consequence of failure is between
114 7.3 Numerical Investigations
Figure 7.6: Likelihood of failure [-] calculated as nf /N versus distance rleaky well between injection
well and leaky well and time t.
zero (no leakage occurred in any case) and 2.72·107 kg, where all cases produced considerable
damage after 7000 days at rleaky well = 100 m. The consequence of failure increases with time
and decreases with distance rleaky well . For any given distance rleaky well the consequence of
failure increases quickly once damage has started to occur. For example occurs a consequence
of failure larger than zero for a leaky well in 200 m distance after 91 days of injection. It
increases then to 106 kg at 792 days, and to 107 kg after 5401 days. Another example is a
leaky well in 2000 m distance, a consequence of failure larger than zero occurs after 2409 days
of injection. It increases then to 104 kg at 2806 days, and to 105 kg after 4819 days.
When looking at the product of both terms, risk as shown in Figure 7.8 is the result. Range
of risk is identical to range of the consequence of failure. This is because the consequence
of failure is multiplied with the likelihood of failure, which is between zero and one. Risk
increases with time and decreases with distance rleaky well , as one would expect.
The risk contour line (red line in Figure 7.8) identifies risk equal to 1 g fitted by a power-
function as given in Equation 7.7.
where r is leaky well distance and A and B are power-fitted coefficients given in Table C.1.
The time contour lines (black lines in Figure 7.8) can be used to calculate risk dependent on
leaky well radii r. They are fitted by a exponential function to the risk surface as given in
7 Risk Analysis 115
PN Di
Figure 7.7: Consequence of failure [log kg] calculated as i=1 nf versus distance rleaky well between
injection well and leaky well and time t.
Equation 7.8.
9
X
Log(Risk/1kg) = Ai · r i , (7.8)
i=0
where r is leaky well distance and Ai are coefficients given in Table C.2.
In Figure 7.9 risk is shown for few selected leaky well distances and for various primary
parameters. In general, risk for all cases increases with time of injection and with smaller
leaky well radii. For a leaky well radius of 200 m, risk is basically identical in all cases. This
result is expected given the short time and spatial scale, because the results are independent
of the multiphase flow and transport regime in the reservoir and consequently independent
of reservoir parameters. The effect and sensitivity of primary parameters is discussed in the
following, where results are shown for clustered groups of cases having the same parameters
(e.g. all cases with porosity of 0.106, 0.190, and 0.320 in Figure 7.9 (top-left)).
116 7.3 Numerical Investigations
PN
Figure 7.8: Log (Risk/1 kg) calculated as nNf · Di
i=1 nf versus distance rleaky well between injection
well and leaky well and time t.
7 Risk Analysis 117
8 All cases with Porosity [-] 8 All cases with Geothermal Gradient [°C]
0.106 0.0204
0.190 0.0323
0.320 0.0532
7 7
200m 200m
6 6
1000m
Log (Risk / 1kg)
4 4
3 3
2000m
2 2
1 1
0 0
0 1000 2000 3000 4000 5000 6000 7000 0 1000 2000 3000 4000 5000 6000 7000
t [d] t [d]
8 All cases with Depth [m] 8 All cases with Anisotropy [-]
640 0.00898
1403 0.01282
2943 0.04274
7 7
200m 200m
6 6
1000m 1000m
Log (Risk / 1kg)
5 5
2000m
4 4
3 3
2000m
2 2
1 1
0 0
0 1000 2000 3000 4000 5000 6000 7000 0 1000 2000 3000 4000 5000 6000 7000
t [d] t [d]
Figure 7.9: Log (Risk/1kg) versus time [d] for selected rleaky well distances (shown in the boxes).
Shown are clustered groups of primary parameters including porosity (top-left),
geothermal gradient (top-right), depth (bottom-left), and absolute permeability
anisotropy (bottom-right).
118 7.3 Numerical Investigations
Figure 7.9 (top-left) shows risk as a function of porosity. No clear trend of risk as a function
of porosity can be seen. Cases with a porosity of 0.106 produce almost the same risk as the
cases do with a porosity of 0.19, and 0.32 respectively. With increasing porosity one would
expect a retardation of the plume evolution velocity, simply due to increased pore space,
and consequently a later increase of risk, i.e. lower risk at any time. But the permeability
also increases with porosity, since it is coupled by a Kozeny-Carman functionality (cf. Sec-
tion 7.2.2). The increased permeability leads to stronger gravity segregation and hence to
increased plume evolution velocity below the caprock. Note that vertical permeability is cou-
pled to horizontal permeability via the Anisotropy parameter. The effects of the retarded
plume evolution velocity due to increased porosity and stronger gravity segregation due to
increased horizontal and vertical permeability cancel each other out, i.e. risk is independent
of porosity. Small differences in risk are due to the other (random) primary parameters
(e.g. depth).
Figure 7.9 (top-right) shows risk for a given geothermal gradient. It can be seen that risk
increases earlier for higher geothermal gradients. This is because of the lower CO2 density for
a higher reservoir temperature. Gravity segregation is stronger at higher temperatures and
fast lateral plume evolution occurs below the caprock. For rleaky well = 1000 m, risk increases
after 884 days for a thermal gradient of 0.0532 ◦ C/m, after 2287 days for a thermal gradient
of 0.0323 ◦ C/m, and after 3538 days for a thermal gradient of 0.0204 ◦ C/m. Risk is also quite
different after 7000 days, i.e. 4.31 · 105 kg, 8.96 · 105 kg, and 1.2 · 106 kg, corresponding to a
factor of 3.
Figure 7.9 (bottom-left) shows risk as a function of depth below ground surface and shows
that risk increases earlier for shallower depth. This is also due to CO2 density, because at
lower density gravity segregation is stronger leading to fast lateral plume evolution below
the caprock. For rleaky well = 1000 m, risk increases after 884 days for a depth of 640 m, after
1433 days for a depth of 1403 m, and after 2684 days for a depth of 2943 m. However, after
7000 days, risk can be viewed as being identical for all depths.
Figure 7.9 (bottom-right) shows risk as a function of permeability anisotropy. Risk increases
earlier for larger anisotropy factors. The gravity segregation effect is stronger for a large
anisotropy factor. Parameter sensitivity, however, is lower, compared to the sensitivities of
geothermal gradient and depth. For rleaky well = 1000 m, risk increases after 1952 days for an
anisotropy of 0.04274, after 2165 days for an anisotropy of 0.01282, and after 2318 days for
an anisotropy of 0.00898. This difference in time grows with increasing rleaky well .
Summarising, it can be concluded that high risk corresponds to short leaky well distances,
long injection times, high geothermal gradients, high permeability anisotropy, and low depth.
Risk, as it is defined here, does not depend on porosity.
7 Risk Analysis 119
Reservoir geometry +/−: An increased height of the reservoir, a negative dip angle towards
the leaky well, and a sealing fracture in the reservoir in between the injection and
the leaky well leads to lower or later leakage through the leaky well. For converse
assumptions, i.e. decreased height, positive dip angle and a sealing fracture on the
opposite side of the injection well, this would lead to higher or earlier leakage. Mosthaf
(2007) studied the influence of varying dip angles (amongst others) and described the
influence by dimensionless numbers. Hesse et al. (2008) suggest that an increasing slope
of an aquifer accelerates residual trapping and that lateral migration of the injected
CO2 traps the CO2 relatively quickly as residual saturation. The results were, however,
obtained with simplified 2-D models neglecting (amongst other simplifications) density
and viscosity changes.
Diffuse leakage through caprock +/−: Diffusive leakage into shallower reservoirs could re-
duce leakage at the leaky well. The effect on risk depends of course on the judgement
whether this diffusive leakage is acceptable (e.g. because leaking CO2 is trapped in the
shallower aquifer) or is not acceptable (leads also to damage or risk). However, Li et al.
(2006) showed, that if the sealing pressure of the caprock is not exceeded, leakage of
CO2 by molecular diffusion is negligible during the short-term injection stage. If the
sealing pressure is exceeded though, leakage rates can become (very) large. Gherardi
et al. (2007) shows that when the transport of chemicals primarily occurs by molecular
diffusion in the water-rich phase, CO2 leakage becomes self-limiting, i.e. pores become
clogged after a very short time. In case of a fractured caprock, however, gaseous CO2
penetrates into the caprock and induces some enhancement in porosity and permeabil-
ity, which reduces the sealing efficiency of the caprock.
plume shape and on dissolved, residually trapped, and mobile CO2 mass fractions.
They conclude that with increasing amounts of shale in the reservoir, vertical move-
ment of the plume is restricted and lateral movement aided. An increasing amount of
shale and consequently a restricted vertical movement thus leads to reduced leakage.
Consequently, it can be concluded that the assumption of homogeneity is a conservative
approach on risk.
Relative Permeability +/−: Considering a different shape of the relative permeability re-
lation or different residual saturations (and consequently different end-points) could
lead to either increased or earlier, or to decreased or later leakage. Bennion and Bachu
(2008) measured relative permeability relations for CO2 -brine systems in the Alberta
basin in Canada, and these do influence the plume evolution behaviour (as shown by
Kopp et al. (2009b)) to a similar degree as the reservoir properties considered here
(depth, geothermal gradient, etc.).
Capillary Pressure +/−: The approach employed to scale a measured capillary pressure
relation to actual reservoir conditions at a given pressure, temperature, and (constant)
porosity is quite sophisticated. However, such a scaling assumes the same rock-type in
all reservoirs. As this is not generally the case, a different capillary pressure relation
might lead to an increase or decrease of leakage. Generally, stronger capillary forces
lead to less leakage, since more CO2 is stored by capillary trapping (Ide et al. (2007),
Kopp et al. (2009b)). In a heterogeneous reservoir, the effect of high capillary entry
pressures of low permeable shale structures amplifies the effects discussed earlier with
respect to heterogeneity, i.e. CO2 may not escape a high-permeability channel structure
towards the leaky well, which leads to earlier leakage and to higher rates.
Anisotropy model +/−: A different anisotropy distribution (here a primary parameter de-
rived by conceptual modeling) could lead to increased or earlier or to a decreased or
later leakage. Generally, an increase in anisotropy leads to decreased or later leakage
due to reduced gravity segregation.
Salinity +/−: As salinity influences brine density, viscosity and the solubility of CO2 in
brine, a higher or lower salinity than assumed influences leakage. Generally, lower
salinity leads to an increase in CO2 solubility in brine (Duan and Sun, 2003), hence
leakage is decreased or occurs later.
Leakage simulation −: If leakage was actually simulated assuming low permeability well
plugs etc. it would always be less than assumed here where leakage composes all the
7 Risk Analysis 121
CO2 mass that passed by the leaky well. In practise, not all of the mass passing by
the leakage point will flow out of the reservoir. Ebigbo et al. (2007) investigated CO2
leakage rates through an abandoned well. The considered reservoir has an open leaky
well (without well plugs etc.) at 100 m distance to the injection well which connects
to a shallower aquifer. An additional simulation is conducted here, resembling the
reservoir parameters and the setup presented by Ebigbo et al. (2007). Hence, the
leakage rate of this additional simulation case derived by the methodology presented
here, and the leakage rate given by Ebigbo et al. (2007) can be compared, as shown in
Figure 7.10.
0.06
0.05
Leakage Rate [%]
0.04
This study
Ebigbo et al. [2007]
0.03
0.0005
Leakage Rate [%]
0.0004
0.02 0.0003
0.0002
0.0001
0.01
0
0 20 40 60 80 100
t [d]
0
0 500 1000 1500 2000
t [d]
Figure 7.10: Comparison of leakage rates versus time obtained by Ebigbo et al. (2007) and cal-
culated by the methodology presented in this study. Leakage rates are given as a
fraction of the CO2 rate injected.
The times when significant leakage commences are almost identical. The difference
is due to slightly different boundary conditions and the different reservoir geometry
(single radially symmetric reservoir in this study versus a fully 3-D simulation with
two reservoirs, connected by the leaky well in the centre). Leakage rates then start to
increase at a slower rate in the study of Ebigbo et al. (2007), reach a peak leakage rate,
and decrease on the long term. The gradual increase of the leakage rate is explained
by up-coning of the brine into the well (Nordbotten et al., 2005b). The higher peak
rate and later decrease is explained by thermodynamic and hydraulic processes in the
leaky well. Pruess (2008a) investigated physical and chemical processes in CO2 leakage
through a fault zone into a shallow reservoir, and eventually consecutive leakage to the
atmosphere along another fault. He outlined some potential for self-enhancement of
122 7.5 Summary and Conclusion
leakage rates by factors of 2–3 relative to CO2 flow rates at depth. This corresponds to
three-phase conditions (water-rich phase – liquid CO2 -rich phase – gaseous CO2 -rich
phase) at relatively shallow depth. But since (i) this happens only temporarily, (ii)
the system was purposefully designed to facilitate such behaviour, and (iii) CO2 was
injected right into the fault, neglecting the processes occurring in the reservoir, the
self-enhancement of leakage is not of importance here. Assuming any kind of (untight)
well plugs etc., reduces leakage rates largely, so that the approach in this study can be
considered a worst case.
Additional assumptions on leaky well properties −: Here the leaky well is assumed to
penetrate the entire formation thickness and is assumed to be completely open to the
atmosphere. In practise this is not the case. For example, partial penetration of the
formation, well plugs, incorporation of detailed leakage pathways in the leaky well,
knowledge about material behaviour etc. will always reduce leakage rates.
CO2 injection scheme and rate +/−: In this study a continuous injection of CO2 is as-
sumed. By a modified injection scheme, i.e. alternating injection of gas and fresh
brine, additional CO2 could be trapped residually (Ide et al., 2007). Kopp et al.
(2009b) showed that the injection rate also has an considerable influence on the shape
of the plume and therefore on the potential leakage. High injection rates lead to rather
cylindrical plume shapes, resulting in less lateral extent of the plume.
Injection well design −: By intelligent design of the injection well (horizontal setup, screen-
ing in lower regions, etc) leakage could be reduced or retarded.
leaky wells in the surrounding of the injection well can be easily incorporated by summing up
the individual risk for each well. The additional assumption herein is, that the leaky wells do
not influence each other, i.e. no leaky well is located in the radial share of a leaky well that
is closer to the injection well (cf. Figure 7.1). Since other than primary reservoir properties
are fixed or have certain dependencies on primary parameters, different assumptions made
for these parameters or for the reservoir could lead to lower/higher or earlier/later leakage
at the leaky well, and hence to lower/larger risk. This is discussed in Section 7.4.
Given the assumptions, this study adopts a conservative approach concerning the leakage
of CO2 through abandoned wells with regard to most reservoir parameters. However, risk
is underestimated if reservoir geometry is very different (large dip, small height) or if the
relative permeability relations show a high CO2 permeability together with large residual
water-rich phase saturation (significantly larger than 0.3).
Main findings and conclusions are :
• The presented study defines a quantitative method for the evaluation of risk with
respect to CO2 leakage through pre-existing wells based on comprehensive reservoir
properties statistics.
• Within the given framework, a range of possible risks is defined. This can be used
to determine whether an individual site is good or not in terms of potential leakage
rate through a number of leaky wells with given distribution around the injection well.
Hence, assistance is given to experts when rating storage sites with unknown/uncertain
reservoir properties. The cumulative risk for any site with given leaky well distances
can be calculated, and sites can be compared to each other. Thus, the relative risk
might assist in making decisions on where to conduct further investigations or to help
experts when utilising more comprehensive screening and ranking frameworks.
• The placement of the injection well can be optimised with respect to risk arising from
abandoned wells in the surroundings. For a given injection well location, the combined
risk for any number of leaky wells in the surroundings can be calculated for the time
of interest analytically. Thus, it is possible to compare risk for several injection well
locations and pre-select the location yielding the lowest risk.
8 Final Remarks
8.1 Summary
Due to tremendous human fossil-fuel use in the past 160 years, the concentration of green-
house gases in the atmosphere increased and is most likely the cause of an observed global
increase of average temperature and of changing climate. It is expected that, with further
global warming, there will be drastic ecological and economic impacts. No one single option
or technology will be sufficient achieve the necessary emission reduction in a growing global
economy.
Carbon dioxide capture and storage considers capturing CO2 in large local sources, such as
power plants etc., transportation to the storage site, and storing the CO2 away from the
atmosphere for geological periods of time. This study focuses on the latter part and con-
siders storage in deep saline geological formations. Injected CO2 spreads in the formation
and is subject to strong buoyancy forces, driving it towards the top of the reservoir where
it is usually kept from further ascent by a low permeable caprock. However, the caprock
may show geological weaknesses or may be perforated by wells which are a potential leakage
pathways for the injected CO2 out of the formation causing risk to health, safety, and the
environment. The careful selection, analysis, and evaluation of a storage site is therefore an
indispensable part of a storage project.
The main objective of this study is the improvement of insight into CO2 injection processes
in geological formations to assist site screening. Site screening is defined as the very first
steps in site selection, where usually little information is available about reservoir properties.
Questions of interest in site screening include the estimation of the storage capacity, which
should be sufficient to store the long-term production of the CO2 source, and the long-term
ability to store CO2 , which is related to the efficiency of the project and risk. These ques-
tions are investigated by the statistical analysis of a database, analytical, and numerical
experiments.
In Chapter 2 the essential features and principal processes of CO2 storage in geological for-
mations are described and a conceptual model is set up. Based on the conceptual model, a
mathematical and numerical model is formulated.
The properties of potential storage sites in geological formations are derived in Chapter 3 by
125
126 8.1 Summary
analysing a large database listing the properties of more than 1200 reservoirs. The parameter
ranges and distributions are used to define typical reservoirs, e.g. a warm, a cold, a deep,
and a shallow reservoir. Additional reservoirs are defined by re-combination of the median
reservoir properties with measured relative permeability-saturations relations, modified cap-
illary entry pressure, halved injection rate and reduced absolute permeability.
With dimensional analysis the dominant forces and relevant processes during the CO2 in-
jection stage are assessed in Chapter 4. This is done by non-dimensionalising the governing
flow equations in the fractional flow formulation by introducing characteristic values for
length, time, velocity, and pressure. The resulting set of balance equations consists only
of dimensionless gradients, dimensionless numbers, and dimensionless saturation-dependent
functions. The dimensionless numbers represent relations of forces in the system, i.e. capil-
lary, viscous, and buoyancy forces. The individual dimensionless terms are then investigated
by independent variation of the characteristic values. This yields a better mathematical un-
derstanding of the system behaviour. To develop a physical understanding of the system, the
mutual dependence of characteristic values and the simultaneous variation of saturations,
gradients, and ratios of forces is investigated in numerical 1-D and 3-D experiments. The
typical reservoirs defined earlier are chosen as a basis.
In Chapter 5, a methodology for investigating the CO2 storage capacity of geological forma-
tions during the injection stage is developed. This approach respects the physical trapping
mechanisms. The storage capacity coefficient C and the effective mass that can be stored
in a reservoir Meff are introduced to compare the storage capacity of the typical reservoirs
defined earlier. To estimate storage capacity, it is necessary to limit the reservoir volume at
some point. This is done by defining a spill point at one kilometre distance from the injection
well. Time-dependent storage capacity estimates are calculated in 1-D and 3-D numerical
experiments. Finally, the results are interpreted using the simultaneously calculated ratios
of forces (dimensionless numbers) in the experiments.
In Chapter 6, the influence of individual reservoir parameters on the model results is inves-
tigated in a sensitivity analysis. The reservoir parameters are varied over their full range
assuming uniform distribution. Model results of interest include measures related to storage
capacity (e.g. arrival time at a spill point) and risk (e.g. developing over-pressure). The ro-
bust, reliable, and efficient Morris Method is used to conduct the sensitivity analysis because
the method is especially suited to screening a large number of model input parameters by
randomly sampling the input parameter space. Qualitative measures are derived identifying
irrelevant input parameters that result in an average input parameter ranking.
In Chapter 7, a risk analysis for potential CO2 leakage through pre-existing leaky wells is per-
formed. Risk is expressed as the likelihood of leakage occurring times the damage produced,
which is defined here as the potentially leaking CO2 mass at a point in time for arbitrary
8 Final Remarks 127
distances between the leaky well and the injection well. The likelihood of leakage occurring
is evaluated by running numerous numerical experiments. In each individual experiment,
properties are randomly sampled from statistical distributions derived from the database
analysis mentioned earlier. The independent parameters are porosity, depth of the reservoir
below surface, average geothermal gradient, and the absolute permeability anisotropy. Other
reservoir parameters are derived by a functional dependency from independent parameters,
or are chosen as constants. In each individual numerical experiment, leakage is calculated
as the accumulated CO2 mass that has passed by the leaky well at a given distance to the
injection well at a given time. Thus, leakage is here primarily subject to multiphase trans-
port processes occurring in the reservoir. The risk surface derived, i.e. risk as a function
of leaky well distance and time, represents average risk for any site with unknown reservoir
properties. The result is discussed comprehensively with respect to the assumptions made.
8.2 Conclusions
The major conclusions and insights into CO2 injection processes in geological formations are
briefly reviewed.
• It is shown that the distributions of analysed reservoir properties do not follow standard
probability distributions.
• The correlation coefficients between investigated properties are low. Mutual parameter
interrelations are rejected, except for the interrelation between absolute permeability
and porosity; here, a Kozeny-Carman relation is assumed.
Dimensional Analysis: Dominant forces and relevant processes during the CO2 injection
stage are assessed by dimensional analysis.
• The dimensionless numbers Ca and Gr are basically able to characterise the plume
evolution during the initial injection phase dominated by multiphase processes.
• The Capillary Number Ca generally reflects plume evolution velocities in 1-D numerical
experiments. Cases with varying injection rate, capillary pressure and permeability do
not fit into this scheme since gravity segregation is neglected as a result of the 1-D
assumption.
not fully understood how the simultaneous variation of Gr and Ca affects estimates of
storage capacity and risk. However, a low ratio of gravitational to viscous forces (low
Gr) - and to some extent also a high ratio of capillary to viscous forces (high Ca)-
possibly leads to high CO2 storage capacity and low risk.
• The relative permeability-saturation relations are of great influence for plume evolution
velocity and average CO2 saturation. Therefore, they have a great influence on storage
capacity and risk.
• The average CO2 saturation can be estimated by analysing the fractional flow function.
• 3-D simulations considering gravity forces show a larger difference in maximum and
minimum storage capacity estimates, than the 1-D simulations.
• In the 3-D simulations, the reservoir with median properties and reduced permeability,
the reservoir with median properties and basal formation relative permeability-relations
and the cold reservoir (in this order) show the highest storage capacity. The warm and
the shallow reservoir shows the smallest storage capacity.
• The estimated storage capacity coefficient C in the 3-D simulations ranges from from
0.0117 (for the shallow reservoir) to 0.036 (for the reservoir with median properties
and reduced permeability). By inclusion of the CO2 density for these setups this
translates to an effective mass that can be stored Meff in a range between 0.09 and
1.94 (normalised to the reservoir with median properties).
• In all 3-D simulations, less than 3.6 % of the total geometric reservoir volume (equiva-
lent to less than 18 % of the pore volume) can be used by the injected CO2 before the
spill point at one kilometre distance is reached.
• This study indicates that kr,CO2 together with the residual water-rich phase saturation
have a major influence on the intrinsic storage capacity coefficient Ci,CO2 , whereas kr,w
has a major influence on the geometric storage capacity coefficient Cg .
8 Final Remarks 129
• Parameters like the salinity, residual water-rich phase saturation, entry pressure, poros-
ity, sorting factor, CO2 injection temperature and the dip angle show the smallest
sensitivity with respect to the questions of interest.
• It is expected that for other setups, like in 3-D simulations, this ordering might be
different due to different relations of the dominant forces in the course of time.
Risk analysis: The risk associated with CO2 leakage through pre-existing wells is analysed
based on comprehensive reservoir-property statistics and a quantitative average risk surface
as a function of the distance of the leaky well distance to the injection well and of time has
been calculated.
• Important reservoir properties with respect to risk are identified based on CO2 plume
evolution behaviour. The reservoir depth below the surface and the geothermal gra-
dient show the highest influence on average risk estimations. Permeability anisotropy
influences the risk estimates only up to some distance from the injection well.
• A large reservoir depth below the surface, a low geothermal gradient, and a large
permeability anisotropy are advantageous for low risk estimates.
• The quantitative average risk surface calculated allows the comparison of individual
sites with respect to risk imposed by the potential leakage through a number of leaky
wells with a given distribution around the injection well. This risk estimate can be
used to decide where to conduct further investigations or help experts utilising more
comprehensive screening and ranking frameworks.
In conclusion, a significant contribution has been made towards the evaluation and the
understanding of CO2 injection processes in geological formations with respect to storage
capacity estimates and risk for leakage through pre-existing wells.
130 8.3 Outlook
8.3 Outlook
In this relatively young, but complex field of research involving various scientific disciplines
and often subject to highly uncertain (geological) data, additional aspects can be included
with the aim of a comprehensive site screening methodology. However, the more aspects are
included in site screening, the more data are necessary which might not be available.
• The aspect of reservoir geometry and property distribution could be included in a di-
mensionless number. A anticlinal shape of the reservoir is advantageous due to limited
spreading of the CO2 . Heterogeneity in the permeability can also be advantageous to
plume spreading, if e.g. gravity segregation is increased. On the other hand, highly
permeable flow paths towards a potential spill point are disadvantageous. The use of
hypothetical sites or real-field project sites would be beneficial here. However, with the
current computational power, it is hardly possible to conduct numerous simulations
(∼100s) of complex reservoir setups for research purposes. With access to computer
clusters or grid environments, this might be different.
• The aspect of an economical CO2 injection into saline aquifers (high injectivity, lower
drilling costs at shallow reservoir depth) could be included by adding another dimen-
sionless number. The injectivity aspect relates to the developing over-pressure, but
other economic aspects may produce a different ranking of the reservoirs here.
• This economic aspect could be further extended by including the aspect of the hy-
drocarbon potential of a site (by enhanced oil recovery, enhanced gas recovery, or
enhanced coal bed methane). However, it would then also be necessary to include all
costs (drilling costs, infrastructure, etc.) when comparing sites.
• More risk-related aspects could be included, e.g. brine displacement to shallower for-
mations, tectonism, etc.
Most of these aspects could be investigated by numerical simulation. Finally, all aspects
could be combined in a comprehensive site screening methodology, allowing for a prelimi-
nary, but quantitative comparison between sites.
Above all, there is a major need for data to develop, apply, and verify techniques and
tools. For this study, more data on measured relative permeability- and capillary pressure-
saturation relations (possibly including hysteresis) would yield great improvements, along
with a global saline-aquifer-properties database to make detailed region-specific storage ca-
pacity estimates possible and to apply and verify developed techniques.
Bibliography
Arrhenius, S., 1896. On the Influence of Carbonic Acid in the Air upon the Temperature of
the Ground. Philosophical Magazine and Journal of Science (5th series) 41, 237–275.
Assteerawatt, A., Bastian, P., Bielinski, A., Breiting, T., Class, H., Ebigbo, A., Eichel,
H., Freiboth, S., Helmig, R., Kopp, A., Niessner, J., Ochs, S., Papafotiou, A., Paul, M.,
Sheta, H., Werner, D., Ölmann, U., 2005. MUFTE-UG: Structure, Applications and Nu-
merical Methods. Newsletter 23 (2), Colorado School of Mine, International Groundwater
Modeling Centre, Golden, Colorado, U.S.A., http://igwmc.mines.edu/news/.
Aziz, K., Settari, A., 1979. Petroleum Reservoir Simulation. Elsevier, Amsterdam, The
Netherlands.
Bachu, S., 2003. Screening and ranking of sedimentary basins for sequestration of CO2 in
geological media in response to climate change. Environmental Geology 44 (3), 277–289.
Bachu, S., Bonijoly, D., Bradshaw, J., Burruss, R., Holloway, S., Christensen, N., Mathi-
assen, O., 2007. CO2 storage capacity estimation: methodology and gaps. International
Journal of Greenhouse Gas Control 1 (4), 430–443.
Baehr, H., Stephan, K., 1998. Wärme und Stoffübertragung. Springer, Berlin, Germany.
Barenblatt, G., Entov, V. M., Ryzhik, V. M., 1990. Theory of Fluid Flows Through Natural
Rocks. Springer, Dordrecht, The Netherlands.
Bastian, P., Birken, K., Johannsen, K., Lang, S., Eckstein, K., Neuss, N., Rentz-Reichert,
H., Wieners, C., 1997. UG - A Flexible Software Toolbox for Solving Partial Differential
Equations. Computing and Visualization in Science, 1(1), 27–40.
Bastian, P., Helmig, R., 1999. Efficient Fully-Coupled Solution Techniques for Two Phase
Flow in Porous Media. Parallel Multigrid Solution and Large Scale Computations. Ad-
vances in Water Resources 23, 199–216.
Batzle, M., Wang, Z., 1992. Seismic Properties of Pore Fluids. Geophysics 57, 1396–1408.
Bear, J., 1972. Dynamics of Fluids in Porous Media. Dover, New York, U.S.A.
Bennion, B., Bachu, S., 2005. Relative Permeability Characteristics for Supercritical CO2
Displacing Water in a Variety of Potential Sequestration Zones in the Western Canada
Sedimentary Basin. Society of Petroleum Engineers, SPE 95547, http://www.spe.org.
133
134 Bibliography
Bennion, B., Bachu, S., 2006. The Impact of Interfacial Tension and Pore Size Distribu-
tion/Capillary Pressure Character on CO2 Relative Permeability at Reservoir Conditions
in CO2 -Brine Systems. Society of Petroleum Engineers, SPE 99325, http://www.spe.org.
Bennion, B., Bachu, S., 2008. Drainage and Imbibition Relative Permeability Relationships
for Supercritical CO2 /Brine and H2 S/Brine Systems in Intergranular Sandstone, Car-
bonate, Shale, and Anhydrite Rocks. SPE Reservoir Evaluation and Engineering 11 (3),
487–496.
Benson, S., 2005. Lessons learned from industrial and natural analogs for health, safety
and environmental risk assessment for geologic storage of carbon dioxide. Carbon Dioxide
Capture for Storage in Deep Geologic Formations - Results from the CO2 Capture Project,
Chapter 25, 1133–1141, Lawrence Berkeley National Laboratory, Berkeley, California,
U.S.A.
Benson, S., Hepple, R., Apps, J., Tsang, C., Lippmann, M., 2005. Lessons Learned from Nat-
ural and Industrial Analogues for Storage of Carbon Dioxide in Deep Geological Forma-
tions. Energy Citations Database, Report Number LBNL–51170, Lawrence Berkeley Na-
tional Laboratory, Berkeley, California, U.S.A., http://repositories.cdlib.org/lbnl/LBNL-
51170/.
Brooks, A., Corey, A., 1964. Hydraulic Properties of Porous Media. Colorado State Univer-
sity Hydrology Paper No.3, Fort Collins, Colorado, U.S.A.
Buckingham, E., 1914. On Physically Similar Systems: Illustrations of the Use of Dimen-
sional Analysis. Physical Review 4, 107–116.
Buckley, S., Leverett, M., 1942. Mechanism of Fluid Displacements in Sands. Transactions
of the AIME 146, 107–116.
Burdine, N., 1953. Relative Permeability Calculations from Pore–Size Distribution Data.
American Institute of Mining, Metallurgical and Petroleum Engineers, Petroleum Trans-
actions 198, 71–77.
Campolongo, F., Braddock, R., 1999. The use of graph theory in the sensitivity analysis of
the model output: a second order screening method. Reliability Engineering and System
Safety 64, 1–12.
Campolongo, F., Cariboni, J., Saltelli, A., Schouters, W., 2004. Enhancing
the Morris Method. Hanson, M.K. and Hemez, F.M. (Editors): 4th Interna-
tional Conference on Sensitivity Analysis of Model Output, Los Alamos National
Bibliography 135
Campolongo, F., Tarantola, S., Saltelli, A., 1999. Tackling quantitatively large dimensional-
ity problems. Computer Physics Communications 117, 75–85.
Celia, M., Bachu, S., 2003. Geological Sequestration of CO2 : Is Leakage Unavoidable and
Acceptable? Gale, J. and Kaya, Y. (Editors): 6th International Conference on Greenhouse
Gas Control Technologies (GHGT-6), Kyoto, Japan, Sept. 30st –Oct. 4th , 477–482.
Celia, M., Bachu, S., Nordbotten, J., Gasda, S., Dahle, H., 2004. Quantitative Estimation of
CO2 Leakage from Geological Storage: Analytical Models, Numerical Models, and Data
Needs. Rubin, E.S. and Keith, D.W. and Gilboy, C.F. (Editors): 7th International Con-
ference on Greenhouse Gas Control Technologies (GHGT-7), Vancouver, Canada, Sept.
5th –9th , p. 9, http://www.ieagreen.org.uk/ghgt7.html.
Chadwick, A., Arts, R., Bernstone, C., May, F., Thibeau, S., Zweigel, P., 2008. Best Practice
for the Storage of CO2 in Saline Aquifers - Observations and Guidelines from the SACS and
CO2STORE projects. British Geological Survey Occasional Publication, 14, Nottingham,
U.K., http://www.bgs.ac.uk/.
Chen, Z., Ewing, R., 1997. Comparison of Various Formulatios of Three-Phase Flow in
Porous Media. Journal of Computational Physics 132, 362–373.
Class, H., 2008. Models for Non-Isothermal Compositional Gas-Liquid Flow and Transport
in Porous Media. Habilitation, Institut für Wasserbau, Universität Stuttgart, Germany,
http://elib.uni-stuttgart.de/opus/frontdoor.php?source opus=3847&la=de.
Class, H., Helmig, R., Bastian, P., 2002. Numerical simulation of non–isothermal multiphase
multicomponent processes in porous media. - 1. An efficient solution technique. Advances
in Water Resources 25, 533–550.
Clauser, C., 1992. Permeability of crystalline rocks. EOS Transactions 73, 233.
Clauser, C., 2003. Numerical Simulation of Reactive Flow in Hot Aquifers, SHEMAT and
Processing SHEMAT. Springer, Berlin.
Clauser, C., 2006. Geothermal Energy, In: K. Heinloth (Editor), GroupVIII: Advanced
Materials and Technologies, Vol. 3: Energy Technologies, Subvol. C: Renewable energies
Edition. Landolt-Börnstein, Springer, Heidelberg.
136 Bibliography
Clauser, C., Huenges, E., 1995. Thermal Conductivity of Rocks and Minerals, In: T.J.
Ahrens (Editor), Rock Physics and Phase Relations - a Handbook of Physical Constants,
AGU Reference Shelf Edition. Vol. 3. American Geophysical Union, Washington.
CO2SINK project consortium, 2009. CO2SINK - In-situ R&D laboratory for geological stor-
age of CO2 ,
http://www.co2sink.org.
Cropp, R., Braddock, R., 2002. The New Morris Method: an efficient second-order screening
method. Reliability Engineering and System Safety 78, 77–83.
Darcy, H., 1856. Les Fontaines Publiques de la Ville de Dijon. Dalmont, Paris, France.
Daubert, T., Danner, R., 1989. Physical and thermodynamic properties of pure chemicals:
data compilation. Hemisphere, New York, U.S.A.
Doughty, C., 2007. Modeling geologic storage of carbon dioxide: Comparison of non-
hysteretic and hysteretic characteristic curves. Energy Conversion and Management 48(6),
1768–1781.
Doughty, C., Pruess, K., Benson, S., Hovorka, S., Knox, P., Green, C., 2001. Capacity
investigation of brine-bearing sands of the Frio-Formation for geological sequestration of
CO2 . Paper LBNL-48176, Lawrence Berkeley National Laboratory, Berkeley, Ca, U.S.A.,
http://repositories.cdlib.org/lbnl/LBNL-48176.
Duan, Z., Sun, R., 2003. An Improved Model Calculating CO2 Solubility in Pure Water and
Aqueous NaCl Solutions from 273 to 533 K and from 0 to 2000 bar. Chem. Geol. 193,
257–271.
Ebigbo, A., Class, H., Helmig, R., 2007. CO2 leakage through an abandoned well: problem-
oriented benchmarks. Computational Geosciences 11(2), 103–115.
Ellert, M., Grønager, M., Konstantinov, A., Konya, B., Lindemann, J., Livenson, I., Nielsen,
J. L., Niinimäki, M., Smirnova, O., Wäänänen, A., 2007. Advanced Resource Connector
middleware for lightweight computational Grids. Future Generation Computer Systems
23, 219–240.
Fenghour, A., Wakeham, W., Vesovic, V., 1998. The Viscosity of Carbon Dioxide. Journal
of Physical and Chemical Reference Data 27(1), 31–44.
Flett, M., Gurton, R., Taggart, I., 2005. The function of gas-water relative permeability
hysteresis in the sequestration of carbon dioxide in saline formations. Society of Petroleum
Engineers SPE 88485-MS, http://www.spe.org/.
Bibliography 137
Flett, M., Gurton, R., Weir, G., 2007. Heterogeneous saline formations for carbon diox-
ide disposal: Impact of varying heterogeneity on containment and trapping. Journal of
Petroleum Science and Engineering 57(1-2), 106–118.
Friedmann, S., 2007. Site Characterization and Selection Guidelines for Geological Car-
bon Sequestration. Report for U.S. Department of Energy by Lawrence Livermore
National Laboratory, Livermore, California, U.S.A., Contract No. W-7405-Eng-48,
http://www.osti.gov/bridge//product.biblio.jsp?query id=0&page=0&osti id=915602.
Garcia, J., 2001. Density of Aqueous Solutions of CO2 . Tech. rep., LBNL Re-
port 49023, Lawrence Berkeley National Laboratory, Berkeley, California, U.S.A.,
http://repositories.cdlib.org/lbnl/LBNL-49023/.
Garcia, J., Pruess, K., 2003. Flow Instabilities During Injection of CO2 into Saline Aquifers.
Tech. rep., Lawrence Berkeley National Laboratory, California, U.S.A., Paper LBNL-
57973, http://repositories.cdlib.org/lbnl/LBNL-57973.
Garcia, J. E., 2003. Fluid Dynamics of Carbon Dioxide Disposal into Saline Aquifers. Ph.D.
thesis, Lawrence Berkeley National Laboratory, California, U.S.A.
Gasda, S., Bachu, S., Celia, M., 2004. The potential for CO2 leakage from storage sites in ge-
ological media: analysis of well distribution in mature sedimentary basins. Environmental
Geology 46(6–7), 707—-720.
Gaus, I., Audigane, P., André, L., Lions, J., Jacquemet, N., Durst, P., Czernichowski-Lauriol,
I., Azaroual, M., 2008. Geochemical and solute transport modelling for CO2 storage, what
to expect from it? International Journal of Greenhouse Gas Control 2(4), 605–625.
Gherardi, F., Xu, T., Pruess, K., 2007. Numerical modeling of self-limiting and self-enhancing
caprock alteration induced by CO2 storage in a depleted gas reservoir. Chemical Geology
244, 103–129.
Gudmundsson, J., Thrainsson, H., 1989. Power potential of two-phase geothermal wells.
Geothermics 18(3), 357–366.
Gunter, W., Wiwchar, B., Perkins, E., 1997. Aquifer disposal of CO2 -rich greenhouse gases:
extension of the time scale of experiments for CO2 -sequestering reactions by geochemical
modelling. Mineralogy and Petrology 59, 121–140.
Helmig, R., 1997. Multiphase flow and transport processes in the subsurface. Springer, Berlin,
Germany.
Helmig, R., Class, H., Huber, R., Sheta, H., Ewing, R., Hinkelmann, R., Jakobs, H., Bastian,
P., 1998. Architecture of the Modular Program System MUFTE UG for Simulating Mul-
tiphase Flow and Transport Processes in Heterogeneous Porous Media. Mathematische
Geologie 2, 123–131.
138 Bibliography
Hesse, M., Orr, F., Tchelepi, H., 2008. Gravity Currents with Residual Trapping. Journal of
Fluid Mechanics 611, 35–60.
Hilfer, R., Øren, P., 1996. Dimensional analysis of pore scale and field scale immiscible
displacement. Transport in Porous Media 22, 53–72.
Holtz, M., 2002. Residual Gas Saturation to Aquifer Influx: A Calculation Method for
3-D Computer Reservoir Model Construction. SPE 75503, paper presented at the SPE
Gas Technology Symposium held in Calgary, Alberta, Canada, April 30st –May 2nd ,
http://www.spe.org/.
Hovorka, S., Benson, S., Doughty, C., Freifeld, B., Sakurai, S., Daley, T., Kharaka, Y., Holtz,
M., Trautz, R., Seay Nance, H., Myer, L., Knauss, K., 2006. Measuring Permanence of
CO2 Storage in Saline Formations: the Frio Experiment. Environmental Geosciences 13,
105–121.
Huber, R., Helmig, R., 1999. Multiphase Flow in Heterogeneous Porous Media: A Classical
Finite Element Method versus an Implicit Pressure - Explicit Saturation-based Mixed
Finite Element - Finite Volume Approach. International Journal for Numerical Methods
in Fluids 29, 899–920.
IAPWS, 1997. Industrial Formulation 1997 for the Thermodynamic Properties of Wa-
ter and Steam, International Association for the Properties of Water and Steam.
http://www.iapws.org/.
Ide, S., Jessen, K., Orr, F., 2007. Storage of CO2 in saline aquifers: Effects of gravity,
viscous, and capillary forces on amount and timing of trapping. International Journal of
Greenhouse Gas Control 1, 481–491.
International Energy Agency Greenhouse Gas R&D Programme, 2008. Aquifer Storage –
Development Issues. Report No. 2008/12, Cheltenham, U.K.
IPCC, 2005. IPCC Special Report on Carbon Dioxide Capture and Storage. Metz, B., David-
son, O., de Coninck, H.C., Loos, M., Meyer, L.A. (Editors): Prepared by Working Group
III of the Intergovernmental Panel on Climate Change, Cambridge University Press, Cam-
bridge, U.K., http://www.ipcc.ch/ipccreports/special-reports.htm.
Johannsen, K., Oswald, S., Held, R., Kinzelbach, W., 2006. Numerical simulation of three-
dimensional saltwater-freshwater fingering instabilities observed in a porous medium. Ad-
vances in Water Resources 29 (11), 1690–1704.
Johannsen, K., Tourunen, O., Kleist, J., 2008. The CO2-Community Grid.
http://wiki.ndgf.org/display/ndgfwiki/CO2.
Bibliography 139
Joint Research Centre of the European Commission, 2004. SIMLAB package Version 2.2, a
simulation environment for sensitivity analysis. Econometrics and Applied Statistics Unit
(EAS), Ispra, Italy, http://simlab.jrc.ec.europa.eu.
Juanes, R., Spiteri, E., Orr, F., Blunt, M., 2006. Impact of relative permeability hysteresis
on geological CO2 storage. Water Resources Research 42 (12).
Kaplan, S., 1997. The Words of Risk Analysis. Risk Analysis 17 (4), 407–408.
Kobus, H., Barczewski, B., Koschitzky, H., 1996. Groundwater and Subsurface Remediation.
Springer, Berlin.
Kopp, A., Class, H., Helmig, R., 2007. Sensitivity Analysis of CO2 Injection Processes in
Brine Aquifers. Presentation given at European Geosciences Union (EGU) General As-
sembly, April 18th , Vienna, Austria.
Kopp, A., Class, H., Helmig, R., 2009a. Investigations on CO2 storage capacity in saline
aquifers - Part 1: Dimensional analysis of flow processes and reservoir characteristics.
International Journal of Greenhouse Gas Control, 3(3), 263–276.
Kopp, A., Class, H., Helmig, R., 2009b. Investigations on CO2 storage capacity in saline
aquifers - Part 2: Estimation of storage capacity coefficients. International Journal of
Greenhouse Gas Control, 3(3), 277–287.
Kühn, M., Clauser, C., 2006. Mineralische Bindung von CO2 bei der Speicherung im Un-
tergrund in geothermischen Reservoiren. Chemie Ingenieur Technik 78(4), Special Issue:
Kohlendioxid und Klimaschutz, 425–434.
Kvamme, B., Kuznetsova, T., Hebach, A., Oberhof, A., Lunde, E., 2007. Measurements and
modelling of interfacial tension for water+carbon dioxide systems at elevated pressures.
Journal of Computational Materials Science 38, 506–513.
Lake, L., 1989. Enhanced Oil Recovery. Prentice Hall Inc., Old Tappan, New York, U.S.A.
Le Gallo, Y., Trenty, L., Michel, A., Vidal-Gilbert, S., Parra, T., Jeannin, L., 2006. Long-
term flow simulations of CO2 storage in saline aquifers. Gale, J. and Bolland, O. (Editors):
8th International Conference on Greenhouse Gas Control Technologies (GHGT-8), Trond-
heim, Norway, June 19th –22nd , p. 6, http://www.ieagreen.org.uk/ghgt8.html.
LeVeque, R., 1992. Numerical Methods for Conservation Laws. Birkhäuser, Basel, Switzer-
land.
Li, Z., Dong, M., Li, S., Huang, S., 2006. CO2 sequestration in depleted oil and gas reservoirs
— caprock characterization and storage capacity. Energy and Conversion and Managment
47, 1372–1382.
Manthey, S., Hassanizadeh, S., Helmig, R., Hilfer, R., 2008. Dimensional analysis of two-
phase flow including a rate-dependent capillary pressure–saturation relationship. Advances
in Water Resources 31(9), 1137–1150.
Marland, G., Andres, B., Boden, T., 2008. Global CO2 Emissions from Fossil-Fuel Burn-
ing, Cement Manufacture, and Gas Flaring: 1751-2005. The Carbon Dioxide Information
Analysis Center (CDIAC), Oak Ridge National Laboratory, Oak Ridge, Tennessee, U.S.A.,
http://cdiac.ornl.gov.
Maul, P., Metcalf, R., Pearce, J., Savage, D., West, J., 2007. Performance assessment for the
geologic storage of carbon dioxide: Learning from the radioactive waste disposal experi-
ence. International Journal of Greenhouse Gas Control 1, 444–455.
Miller, C., Christakos, G., Imhoff, P., McBride, J., Pedit, J., Trangenstein, J., 1996. Mul-
tiphase flow and transport modeling in heterogeneous porous media: challenges and ap-
proaches. Advances in Water Resources 21(2), 77–120.
Mito, S., Xue, Z., Ohsumi, T., 2008. Case study of geochemical reactions at the Nagaoka
CO2 injection site, Japan. International Journal of Greenhouse Gas Control 2, 309–318.
Morris, M., 1991. Factorial Sampling Plans for Preliminary Computational Experiments.
Technometrics 33(2), 161–174.
Mosthaf, K., 2007. CO2 Storage into Dipped Saline Aquifers Including Ambient Water Flow.
Diploma Thesis, Institute of Hydraulic Engineering, University Stuttgart, Stuttgart, Ger-
many, http://www.iws.uni-stuttgart.de/publikationen/ausgabe.php?person=1&autor[0]=
1195.
Nagaoka project consortium, 2009. Research Institute of Innovative Technogoly for the
Earth, CO2 Sequestration Research Group, Kizugawadai, Kizugawa-shi, Kyoto, Japan,
http://www.rite.or.jp/English/lab/geological/geological.html.
NETL, 2007. Carbon Sequestration Atlas of the United States and Canada. U.S. Department
of Energy - National Energy Technology Laboratory, Pittsburgh, Pennsylvania, U.S.A.,
http://www.netl.doe.gov/technologies/carbon seq/refshelf/atlasII/index.html.
Nordbotten, J., Celia, M., Bachu, S., 2004. Analytical Solutions for Leakage Rates through
Abandoned Wells. Water Resources Research 40(4), W04204.
Nordbotten, J., Celia, M., Bachu, S., 2005a. Injection and Storage of CO2 in Deep Saline
Aquifers: Analytical Solution for CO2 Plume Evolution During Injection. Transport in
Porous Media 58(3), 339–360.
Nordbotten, J., Celia, M., Bachu, S., Dahle, H., 2005b. Semi-Analytical Solution for CO2
Leakage through an Abandoned Well. Environmental Science and Technology 39(2), 602–
611.
NPC, 1984. U.S. National Petroleum Council Public Database. Washington D.C., U.S.A.,
http://www.netl.doe.gov/technologies/oil-gas/Software/database.html.
Ochs, S., 2006. Steam injection into saturated porous media - process anal-
ysis including experimental and numerical investigations -. Ph.D. thesis, Uni-
versität Stuttgart, Stuttgart, Germany, Mitteilungsheft Nr. 159, http://elib.uni-
stuttgart.de/opus/frontdoor.php?source opus=2971&la=de,.
Oldenburg, C., 2007. Screening and ranking framework for geologic CO2 storage site selection
on the basis of health, safety, and environmental risk. Environmental Geology 54 (8), 1687–
1694.
Pacala, S., 2003. Global constraints on reservoir leakage. Gale, J. and Kaya, Y. (Editors):
6th International Conference on Greenhouse Gas Control Technologies (GHGT-6), Kyoto,
Japan, Sept. 30st –Oct. 4th , 267–272.
Pacala, S., Socolow, R., 2004. Stabilization Wedges: Solving the Climate Problem for the
Next 50 Years with Current Technologies. Science 305, 968–972.
Pape, H., Clauser, C., Iffland, J., 1999. Permeability prediction based on fractal pore-space
geometry. Geophysics 64(5), 1447–1460.
Peng, D.-Y., Robinson, D. B., 1976. A New Two-Constant Equation of State. Industrial and
Engineering Chemistry Fundamentals 15, 59–64.
142 Bibliography
Plug, W.-J., Bruining, J., 2007. Capillary pressure for the sand-CO2 -water system under var-
ious pressure conditions. Application to CO2 sequestration. Advances in Water Resources
30 (11), 2339–2353.
Probst, P., 2008. Numerical Simulations of CO2 Injection into Saline Aquifers:
Estimation of Storage Capacity and Arrival Times using Multiple Realiza-
tions of Heterogeneous Permeability Fields. Diploma Thesis, Institute of Hy-
draulic Engineering, University Stuttgart, Stuttgart, Germany, http://www.iws.uni-
stuttgart.de/publikationen/ausgabe.php?diplom=1&betreuer=484.
Pruess, K., 2008a. Leakage of CO2 from geologic storage: Role of secondary accumulation
at shallow depth. International Journal of Greenhouse Gas Control 2, 37–46.
Pruess, K., 2008b. On CO2 fluid flow and heat transfer behavior in the subsurface, following
leakage from a geologic storage reservoir. Environmental Geology 54, 1677–1686.
Pruess, K., Bielinski, A., Ennis-King, J., Fabriol, R., Le Gallo, Y., Garcia, J., Jessen, K.,
Kovscek, T., Law, D.-S., Lichtner, P., Oldenburg, C., Pawar, R., Rutqvist, J., Steefel, C.,
Travis, B., Tsang, C.-F., White, S., Xu, T., 2003. Code Intercomparison Builds Confidence
in Numerical Models for Geologic Disposal of CO2. Gale, J. and Kaya, Y. (Editors): 6th
International Conference on Greenhouse Gas Control Technologies (GHGT-6), Kyoto,
Japan, Sept. 30st –Oct. 4th , 463–470.
Rutqvist, J., Birkholzer, J., Tsang, C.-F., 2008. Coupled reservoir–geomechanical analysis
of the potential for tensile and shear failure associated with CO2 injection in multilayered
reservoir–caprock systems. International Journal of Rock Mechanics and Mining Sciences
45(2), 132–143.
Saltelli, A., Chan, K., Scott, E. M. (Eds.), 2000. Sensitivity Analysis. John Wiley & Sons,
Chichester, U.K.
Scheidegger, A. E., 1960. The physics of flow through porous media. Macmillan Co., New
York, U.S.A.
Somerton, W., El-Shaarani, A., Mobarak, S., 1974. High Temperature Behaviour of
Rocks Associated with Geothermal Type Reservoirs. SPE-4897-MS, paper presented at
the SPE California Regional Meeting, 4th –5th April, San Francisco, California, U.S.A.,
http://www.spe.org/.
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.
Bibliography 143
Spiteri, E., Juanes, R., Blunt, M., Orr, F., 2005. Relative-Permeability Hysteresis: Trapping
Models and Application to Geological CO2 Sequestration. SPE 96448-MS, paper presented
at the Annual Technical Conference and Exhibition, 9th –12th October, Dallas, Texas,
U.S.A., http://www.spe.org/.
Spiteri, E., Juanes, R., Blunt, M., Orr, F., 2008. A New Model of Trapping and Relative
Permeability Hysteresis for All Wettability Characteristics. Society of Petroleum Engineers
Journal 13(3), 277–288.
Stern, N., 2007. The Economics of Climate Change: The Stern Review. Cambridge University
Press, Cambridge, U.K.
Torp, T., Gale, J., 2004. Demonstrating storage of CO2 in geological reservoirs: The Sleipner
and SACS projects. Energy 29 (9–10), 1361–1369.
United Nations Framework Convention on Climate Change, 1997. The Kyoto Proto-
col. The Climate Change Secretariat, Bonn, Germany, http://unfccc.int/kyoto proto-
col/items/2830.php.
Van Der Waals, J. D., 1873. Over de continuiteit van den gas- en vloeistoftoestand. Ph.D.
thesis, Universiteit Leiden, Leiden, The Netherlands.
Van Genuchten, M., 1980. A closed-form equation for predicting the hydraulic conductivity
of unsaturated soils. Soil Science Society of America Journal 44, 892–898.
Xu, T., Apps, J., Pruess, K., Yamamoto, H., 2007. Numerical modeling of injection and
mineral trapping of CO2 with H2 S and SO2 in a sandstone formation. Chemical Geology,
242, 319–346.
Xu, T., Pruess, K., 1998. Coupled modeling of nonisothermal multiphase flow, solute trans-
port and reactive chemistry in porous and fractured media: 1. Model development and
validation. Report Number LBNL–42050, Lawrence Berkeley National Laboratory, Berke-
ley, California, U.S.A., http://repositories.cdlib.org/lbnl/LBNL-42050/.
A Mathematical and Numerical Model
Z Z
dE ∂(% e)
= dΩ + (% e)(v · n) dΓ, (A.1)
dt ∂t
Ω Γ
where Ω is the control volume domain and Γ is the control volume boundary. The vector n
is normal to the boundary Γ of the control volume domain.
∂%
+ ∇ · (% v) = 0. (A.2)
∂t
Conservation of Momentum
The momentum conservation is set up by inserting (mv) for the extensive property E (con-
sequently e = v). When applying Newton’s second law, i.e. the rate of change of momentum
( dE
P
dt
) is equal to the external forces ( F ) exerted upon the system, and the Green-Gaussian
R R
integral rule Γ (% v)(v · n) dΓ = Ω ∇ · (%v · v) dΩ is applied, the momentum conservation
reads in differential form:
P
d(% v) F
+ ∇· (% v · v) = , (A.3)
dt V
145
146 A.1 Mathematical Model for Multi-Phase Processes - the 2p-module
where V is the volume of the domain Ω. Equation A.3 is valid for describing the momentum
conservation of any fluid on the microscale. On the macroscale however, the multi-phase
extension of the empirically derived Darcy law (Equation 2.30) is used to calculate the phase
velocity vα . The advantage is that the phase velocity is explicitly given and can be inserted
into the mass balance equation (Equation A.4).
∂(φ Sα %α )
+ ∇ · (%α vα ) − qα = 0, α ∈ {w, CO2 }. (A.4)
∂t
Simplifying Assumptions
The accurate solution of any multi-phase system in a porous medium would result in a large
system of balance equations. The computational cost of the solution of such a large system
is very high and simplifying assumptions can be made to reduce the cost. It is of importance
to make these assumptions without reducing the accuracy of the model considerably. In the
context of this study, it can be assumed that the solid phase (the porous media) is rigid
(thus ∂φ
∂t
= 0), immobile, and inert.
Local thermodynamic equilibrium is assumed, i.e. local thermal, mechanical, and chemical
equilibrium (cf. Section 2.2.5). However, thermal and chemical equilibrium is not of impor-
tance in this module, since temperature effects are neglected (no energy balance is solved)
and mass transfer between the phases is also neglected (mutual dissolution of components
A Mathematical and Numerical Model 147
is not considered). Thus, only mechanical equilibrium is valid here as the capillary pressure
difference between the phases on the macroscale is considered. This will be different for the
non-isothermal multi-phase multi-component model (Section A.2).
Equation of State
The equation of state for a pure substance is a mathematical formulation describing the
equilibrium relationship between pressure, temperature, and volume. The equation of state
for the system considered here is discussed in Section 2.1.4.
Constitutive Relations
Typical constitutive relations for multi-phase systems include the capillary pressure-saturation
relations and the relative permeability-saturation relations (cf. Section 2.2.4). These rela-
tions are typically of empirical nature; thus, they are most often only approximate and
uncertain.
Auxiliary Conditions
Auxiliary conditions are necessary to close the system of equations and result directly from
the definitions of saturation (Equation 2.6) and of the capillary pressure relation on the
macroscale (Equation 2.23) already known. This means phase saturations have to add up
to one and the phase pressure of one phase can be calculated from the phase pressure of the
respective other phase and the capillary pressure (which is a function of saturation on the
macroscale).
Conservation of Energy
The first law of thermodynamics states that the energy in a closed system is conserved
(Baehr and Stephan, 1998). Mathematically expressed, this means that the change in the
internal energy (dU ) is equal to the amount of energy added by heating (Q) plus the amount
of energy added by doing work on the system (dW ):
dU = dQ + dW. (A.5)
Applying the Reynolds transport theorem (Equation A.1) for internal energy as the extensive
property, i.e. E = U and e = u, yields:
Z Z
dU ∂(% u) dQ dW
= dΩ + (% u)(v · n) dΓ = + . (A.6)
dt ∂t dt dt
Ω Γ
A heat flux over the system boundaries may result from radiation and heat conduction only.
However, when radiation due to small temperature gradients in the subsurface is neglected,
the heat flux term reads as follows:
Z
dQ
= − (qh · n) dΓ. (A.7)
dt
Γ
R R
Applying the Green-Gaussian integral rule Γ (qh · n) dΓ = Ω ∇ · qh dΩ (Helmig, 1997) and
inserting Fourier’s law (Equation 2.38) for the heat conduction qh leads to:
Z
dQ
= ∇ · (λi ∇T ) dΩ. (A.8)
dt
Ω
The change of internal energy due to the amount of energy added by doing work can either
result from dissipative work or from volume-changing work. Due to small flow velocities
A Mathematical and Numerical Model 149
dissipative work can be neglected and volume-changing work is then expressed as (Ochs,
2006):
Z
dW
= − ∇ · (p v) dΩ. (A.9)
dt
Ω
R R
Following application of the Green-Gaussian integral rule Γ (% u)(v·n) dΓ = Ω ∇·(% u v) dΩ,
and insertion of Equations A.8 and A.9 into Equation A.6, energy conservation is expressed
as:
Z Z Z Z
dU ∂(% u)
= dΩ + ∇ · (% u v) dΩ = ∇· λi ∇T dΩ − ∇ · (p v) dΩ. (A.10)
dt ∂t
Ω Ω Ω Ω
p
When the definition of the specific enthalpy h = u + %
(Equation 2.15) is inserted, energy
conservation can be written in differential form as:
∂(% u)
+ ∇ · (% h v) − ∇· λi ∇T = 0. (A.11)
∂t
P
C
∂ α (%α Xα Sα )
φ
| ∂t
{z }
storage
X
C
− ∇· %α Xα k λα (∇pα − %α g ∇z)
|α {z }
advective transport
− ∇· DC C
pm %w ∇Xw
| {z }
diffusive transport
The different terms can clearly be identified, i.e. a storage term, an advective transport
term, a diffusive transport term, and a term considering sources and sinks.
P
∂ (% u S
α α α α ) ∂(%s cs T )
φ + (1 − φ)
| ∂t
{z } | {z ∂t }
energy storage in the fluids energy storage in the matrix
− ∇ · (λpm ∇T )
| {z }
heat conduction
X
− ∇· %α hα k λα (∇pα − %α g ∇z)
|α {z }
heat transport due to advection
X
− ∇· DC % hC
pm w w ∇X C
w
|C {z }
heat transport due to diffusion
where %s is the solid grain density, cs is the specific heat capacity of the soil grains, λpm is the
local heat conductivity as a function of the heat conductivity of the matrix and the fluids, of
porosity, and of the phase saturations (cf. Section 2.3.5), and DC pm is the diffusion coefficient
of the components in the porous medium.
Simplifying Assumptions
In addition to the mechanical equilibrium valid for the 2p-module (considering the capillary
pressure difference between the phases on the macroscale) in the mathematical model for
non-isothermal multi-phase multi-component processes, local thermal and chemical equilib-
rium is of importance (cf. Section 2.2.5). Assuming local thermal equilibrium allows the
solution of just one energy balance for the system, instead of one energy balance for each
phase. Assuming local chemical equilibrium allows the calculation of the mass fraction of a
component in all phases by knowing the mass fraction of a component in one phase (Miller
et al., 1996). Thus, any mass transfer kinetics are neglected and mass fractions of compo-
nents are instantaneously in equilibrium among the phases.
Equation of State
Equations of state cannot only be defined for pure substances (cf. Section A.1.3), but also
for mixtures of substances. To describe mixtures, mixing rules have to be set up, taking
into account the properties of the pure substances and the interaction effects between the
substances. This is discussed in Section 2.1.4.
Constitutive Relations
In addition to the constitutive relations given in Section A.1.3, rules for the mass transfer
between the phases need to be found for the multi-component model considered here. This
is discussed in Section 2.3.7.
Auxiliary Conditions
To close the system of equations for the multi-component model considered here, one addi-
tional auxiliary condition is necessary. The condition follows directly from the definition of
the mass fractions of the components in the phases, i.e. mass fractions have to add up to
one in each phase (Equation 2.3).
Substitution criteria
Phase state Present ph. Primary variables
W-rich phase appears CO2 -rich phase appears
Both phases w, CO2 pCO2 , T , Sw – –
W-rich phase w pCO2 , T , XwCO2 – XwCO2 ≥ (XwCO2 )max
w w w
CO2 -rich ph. CO2 pCO2 , T , XCO2 XCO2 ≥ (XCO2 )max –
Table A.1: Primary variables and substitution criteria for the 2p2cni-module (Bielinski, 2006).
A detailed description of the process-adaptive algorithm for the substitution of the primary
variables and the implementation in the numerical simulator used here is given in Class et al.
(2002) and Bielinski (2006).
Dirichlet type: The Dirichlet boundary condition defines the value of a primary variable at
w
the boundary of the model domain (e.g. Sw , pCO2 , T, XCO2 , XwCO2 ).
Neumann type: The Neumann boundary condition describes the gradient of a primary
variable at the boundary; it is therefore a flux of a quantity perpendicular to the
boundary of the model domain. Here, this is a mass or energy flux.
Integrating the differential form of the mass balance equation A.12 over the model domain
yields the weak form of the equations:
P
C
Z ∂ α (%α Xα Sα )
φ dΩ
∂t
Ω
Z X
− ∇· %α XαC k λα (∇pα − %α g ∇z) dΩ
Ω α
Z
− ∇· DC %
pm w ∇X C
w dΩ
Ω
Z
− q C dΩ = 0, C ∈ {w, CO2 }, α ∈ {w, CO2 }. (A.14)
Ω
w
The discrete values of the primary variables u (e.g. Sw , pCO2 , T, XCO2 , XwCO2 ) are given at the
nodes of the finite element mesh. In between these nodes, the values are approximated using
a basis function Nj for node j:
nX
nodes
ũ = ûj · Nj , (A.15)
j=1
where nnodes is the number of nodes of the finite element mesh. Figure A.1 shows the basis
function, which is a C 0 Lagrangian polynomial.
Ni N i+1
1.0
0.0
i+2
x
i−1 i i+1
Figure A.1: Basis function for the respective node in the 1-D case. The basis function is one at
the respective node and zero at all other nodes.
Inserting the approximated values of the primary variables into the balance equations results
in an error ε. The weighted mean error shall become zero over the model domain.
Z
!
Wi · ε dΩ = 0, i = 1, 2, . . . , nnodes , (A.16)
Ω
154 A.4 Discretisation in Space and Time
Following application of the implicit Euler scheme for the time discretisation and inclusion
of the definitions made in Equations A.15 to A.18, Equation A.14 is then rewritten as
!
φ X X (t+∆t) X t
Mij (%α X̂αC Ŝα ) − (%α X̂αC Ŝα )
∆t j∈η α
j
α
j
i
Z XX
(t+∆t)
− Wi ∇ · (%α k λα )ij (X̂αC Nj )(t+∆t) ∇Nj dΩ (Ψ̂αj − Ψ̂αi )(t+∆t)
Ω α j∈ηi
Z X
(t+∆t)
− Wi ∇· (DC
pm %w )ij C
∇Nj dΩ (X̂wj C (t+∆t)
− X̂wi )
Ω j∈ηi
Z
− (Wi q̂iC )(t+∆t) dΩ = 0, (A.19)
Ω
where ηi is a set of nodes including all neighbouring nodes of the considered node i.
A mass lumping technique is introduced, assigning all entries of the mass matrix to its main
diagonal (Huber and Helmig, 1999):
lumped Vi for i = j
Mij = (A.20)
0 for i 6= j
R
Additionally, the source term is lumped by Ω Wi qiC dΩ = Vi qiC .
R R R
Applying the product rule Ω ∇ · (Wi · F) dΩ = Ω (∇Wi · F) dΩ + Ω (Wi ∇ · F) dΩ and applying
R R
the Green-Gaussian integral rule (Helmig, 1997) as Ω ∇ · (Wi · F) dΩ = Γ (Wi · F) · n dΓ
yields after rearrangement:
Z Z Z
(Wi ∇ · F) dΩ = (Wi · F) · n dΓ − (∇Wi · F) dΩ. (A.21)
Ω Γ Ω
A Mathematical and Numerical Model 155
Equation A.21 is introduced exemplarily to the advective flux term in Equation A.19:
Z XX
(t+∆t)
Wi ∇ · (%α k λα )ij (X̂αC Nj ) (t+∆t)
∇Nj dΩ (Ψ̂αj − Ψ̂αi )(t+∆t)
Ω α j∈ηi
Z XX !
(t+∆t)
= Wi · (%α k λα )ij (X̂αC Nj ) (t+∆t)
∇Nj n dΓ (Ψ̂αj − Ψ̂αi )(t+∆t)
α j∈ηi
ZΓ X X
(t+∆t)
− ∇Wi (%α k λα )ij (X̂αC Nj )(t+∆t) ∇Nj dΩ (Ψ̂αj − Ψ̂αi )(t+∆t) (A.22)
Ω α j∈ηi
For the discretisation in space, the Box method (Helmig (1997), Bastian and Helmig (1999))
is used; this is a node-centred (vertex-centred) finite volume method based on the Galerkin
finite element method. To create the primary finite element mesh, the model domain is split
into elements that connect at the nodes of the mesh (see Figure A.2). A secondary mesh
is constructed by connecting the midpoints of the element edges with the barycentres of
the elements. Thus, each element is split into a number of sub-control volumes. The sub-
control volumes around the respective point form a box. The fluxes between the boxes are
approximated at the integration points, which lie at the centre of each sub-control volume
face.
box B i
o
n
p
sub control volume (SCV)
Bi
m
i i
j j
e
e
l k barycenter of k integration points (IPs)
element e
Figure A.2: Finite element mesh (solid lines) and finite volume mesh (dashed lines) composed of
several sub control volumes (Helmig, 1997).
The Box method employs a weighting function as shown in Figure A.3. Since the weighting
function is constantly one inside the respective box, the gradient ∇Wi equals zero. Hence,
the last term in Equation A.22 vanishes.
To obtain a stable and non-oscillating solution, it is necessary to evaluate the mobilities
λα , the densities %α , and mass fractions of the components XαC in the advective term at
the upstream node (Helmig, 1997). The upstream node can be identified by its larger total
potential, compared to the total potential of the downstream node. This is referred to as
156 A.4 Discretisation in Space and Time
Wi Wi+1
1.0
0.0
x
i−1 i i+1
Figure A.3: Weighting function (W) for the respective node in the 1-D case. For the box method
discussed here, the weighting function is one inside the respective box and zero else-
where.
“fully-upwinding”:
i if (Ψαj − Ψαi ) ≤ 0
upw(i, j) = (A.23)
j if (Ψαj − Ψαi ) > 0
The coefficients in the diffusive flux term are evaluated by arithmetic averaging between
adjacent nodes.
Finally, after inclusion of the mass lumping (Equation A.20), application of the Green-
Gaussian integral rule (Equation A.21) to the advective and diffusive flux terms (shown ex-
emplarily for the advective flux term in Equation A.22), introduction of the fully-upwinding
technique (Equation A.23) and on the assumption weighting functions as shown in Fig-
ure A.3, Equation A.19 is rewritten as:
!
φ X (t+∆t) X t
Vi (%α X̂αC Ŝα ) − (%α X̂αC Ŝα )
∆t α
j
α
j
Z X X
− (%α λα X̂αC )(t+∆t)
upw k ∇Nj n dΓ (Ψ̂αj − Ψ̂αi )(t+∆t)
Γ α j∈ηi
Z X
(t+∆t)
− (DC
pm %w )ij
C
∇Nj n dΓ (X̂wj C (t+∆t)
− X̂wi )
Γ j∈ηi
R(u) = 0, (A.25)
−1 (t+∆t,n)
u(t+∆t,n+1) = u(t+∆t,n) − Jt+∆t,n ·R u . (A.27)
The Jacobian matrix J can be calculated using a central difference scheme as:
For the solution of the linearised equations, MUFTE-UG offers various sophisticated meth-
ods, e.g. the bi-conjugated gradient method (Bastian et al., 1997) or the multi-grid method
(Class et al., 2002).
B Detailed derivation of Equations 4.21,
4.22, and 4.25
Inserting Equations 4.12 to 4.17 into the pressure and water-phase saturation Equations
(Equations 4.8, 4.9 and 4.10) leads to a dimensionless pressure equation
∇·
b v̂tot = 0, (B.1)
1
bpw pcr + fCO2 1 ∇bbpc pcr − g 1 ∇ẑ
X
v̂tot · vcr = − λ k ∇c b lcr fα %α , (B.2)
lcr lcr lcr α
(B.3)
∂Sw 1
φ + ∇
b · fw v̂tot · vcr (B.4)
∂ t̂ · tcr lcr
1 1 b
+ ∇ · fw λCO2 (%w − %CO2 ) k g ∇ẑ lcr
b
lcr lcr
1 1
+ ∇
b · fw λCO2 k ∇b bpc pcr = 0.
lcr lcr
λCO2
Some reformulation (using e.g. fCO2 = λ
) yields for the pressure equation
∇·
b v̂tot = 0, (B.5)
"
bpw pcr λ k + ∇b
v̂tot = − ∇c bpc pcr λCO2 k (B.6)
lcr vcr lcr vcr
#
g k λ lcr
%w f CO2 % CO2 f w
− ∇ẑ
b (%w − %CO2 ) (fw − fCO2 ) + + ,
vcr lcr %w − %CO2 %w − %CO2
| P
{z }
fα % α
α
159
160
tcr b
+ ∇ · fw λCO2 (%w − %CO2 ) k g ∇ẑ
b
φ lcr
|{z}
v−1
cr
tcr b 1 b
+ ∇ · fw λCO2 k ∇b pc pcr = 0.
φ lcr lcr
|{z}
v−1
cr
kr,CO2
Some further reformulation (using e.g. λ = 1
µCO2
( µµCO2
w
kr,w + kr,CO2 ), and λCO2 = µCO2
) yields
for the pressure equation
∇·
b v̂tot = 0, (B.8)
"
k pcr µ
CO2
k pcr
v̂tot = − ∇c
bpw kr,w + kr,CO2 + ∇b
bpc kr,CO2 (B.9)
µCO2 vcr lcr µw µCO2 vcr lcr
| {z } | {z } | {z }
Ca A Ca
Table C.1: Power-fitted coefficients A and B in Equation 7.7 to calculate risk contour lines. R2
indicates goodness of fit (that is the ratio of the sum of the squares of the regression
to the total sum of the squares) between {0,1}, where a value closer to one indicates
better fit.
161
162
Time A0 A1 A2 A3 A4 A5 A6 A7 A8 A9 R2
contour
line
[days]
1000 1.0685E+00 1.4859E-01 -1.5375E-03 8.0130E-06 -2.3983E-08 4.2862E-11 -4.5849E-14 2.8163E-17 -8.8019E-21 9.8655E-25 0.9997
2000 5.0088E+00 4.9038E-02 -4.6774E-04 2.0879E-06 -5.2227E-09 7.7607E-12 -6.9974E-15 3.7556E-18 -1.1029E-21 1.3634E-25 0.9998
3000 7.4770E+00 -3.9704E-03 -7.6240E-06 7.2986E-08 -1.9675E-10 2.5317E-13 -1.7582E-16 6.6951E-20 -1.2949E-23 9.6473E-28 0.9995
4000 8.4316E+00 -1.9807E-02 9.8945E-05 -2.7483E-07 4.3676E-10 -4.2399E-13 2.5517E-16 -9.2601E-20 1.8511E-23 -1.5627E-27 0.9987
5000 8.0576E+00 -1.0807E-02 3.9636E-05 -8.7095E-08 1.1024E-10 -8.6559E-14 4.2855E-17 -1.2960E-20 2.1768E-24 -1.5516E-28 0.9973
6000 8.1624E+00 -1.0943E-02 3.8833E-05 -8.0753E-08 9.6220E-11 -7.0288E-14 3.1944E-17 -8.7698E-21 1.3262E-24 -8.4625E-29 0.9974
7000 8.3008E+00 -1.1898E-02 4.2956E-05 -8.8591E-08 1.0433E-10 -7.4669E-14 3.2931E-17 -8.7074E-21 1.2621E-24 -7.6946E-29 0.9964
Table C.2: Coefficients Ai fitted to Equation 7.8 to calculate time contour lines. R2 indicates goodness of fit (that is the ratio of the
sum of the squares of the regression to the total sum of the squares) between {0,1}, where a value closer to one indicates
better fit.
D Output List
D.1 Peer-Reviewed
• Kopp, A., Binning, P.J., Johannsen, K., Class, H. and Helmig, R., Risk Analysis for
Leakage through Abandoned Wells in Geological CO2 Storage, submitted to Energy
Conversion and Management, 2009.
• Kopp, A., Ebigbo, A., Bielinski, A., Class, H. and R. Helmig, Numerical simulation
of temperature changes due to CO2 injection in geological reservoirs, AAPG Studies
59: Carbon dioxide sequestration in geological media – State of the science, Grobe, M.
and Pashin, J.C. and Dodge, R.L. (Editors), in press, 2009.
• Kopp, A., Class, H. and Helmig, R., Investigations on CO2 storage capacity in saline
aquifers - Part 1: Dimensional analysis of flow processes and reservoir characteristics,
Int. J. Greenhouse Gas Control, 3(3), 263–276, DOI:10.1016/j.ijggc.2008.10.002, 2009.
• Kopp, A., Class, H. and Helmig, R., Investigations on CO2 storage capacity in saline
aquifers - Part 2: Estimation of storage capacity coefficients, Int. J. Greenhouse Gas
Control, 3(3), 277–287, DOI:10.1016/j.ijggc.2008.10.001, 2009.
• Bielinski, A., Kopp, A., Schütt, H. and Class, H., Monitoring of CO2 plumes during
storage in geological formations using temperature signals: numerical investigation,
Int. J. Greenhouse Gas Control, 2, 319–328, DOI: 10.1016/ j.ijggc.2008.02.008, 2008.
• Hurter, S., Garnett, A., Bielinski, A. and Kopp, A.,Thermal Signature of Free-Phase
CO2 in Porous Rocks: Detectability of CO2 by Temperature Logging, Society of
Petroleum Engineers, SPE 109007, DOI: 10.2118/109007-MS, 2007.
• Class, H., Bielinski, A., Helmig, R., Kopp, A. and Ebigbo, A., Numerische Simulation
der Speicherung von CO2 in geologischen Formationen, Chemie Ingenieur Technik,
78(4), 445–452, Special Issue: Kohlendioxid und Klimaschutz, DOI: 10.1002/cite.200500186,
2006.
163
164 D.2 Non Peer-Reviewed (selected)
• Johannsen, K., Kopp, A., Tourunen, O. and Kleist, J., Studying CO2 Sequestration
with the Power of Supercomputing, European Research Consortium for Informatics
and Mathematics (ERCIM) 74, 17–18, 2008.
• Class, H., Ebigbo, A. and Kopp, A., Numerical modeling of CO2 storage in geological
formations - recent developments and challenges, 1. French-German Symposium on
Geological Storage of CO2 , Science Report No.9, Geoforschungszentrum Potsdam, 51–
58, 2007.
• Kopp, A., Class, H. and Helmig, R., Treibhausgase in der Bodenfalle, UNI-Kurier
1/2007,99, EU-Projekt CO2 SINK: Numerische Simulation der Kohlendioxidspeicherung,
2007.
• Kopp, A., Bielinski, A., Ebigbo, A., Class, H. and Helmig, R., Numerical Inves-
tigation of Temperature Effects during the Injection of Carbon Dioxide into Brine
Aquifers, Gale, J. Bolland, O. (Editors): 8th International Conference on Green-
house Gas Control Technologies, (GHGT-8), Trondheim, Norway, June 19th –22nd ,
http://www.ieagreen.org.uk/ghgt8.html.
• Kopp, A. and Sheta, H., Grid Generation for Complex Geological Systems in Mining
Areas, European Union Marie Curie Conferences and Training Courses, 2006.
• Assteerawatt, A., Bastian, P., Bielinski, A., Breiting, T., Class, H., Ebigbo, A., Eichel,
H., Freiboth, S., Helmig, R., Kopp, A., Niessner, J., Ochs, S.O., Papafotiou, A., Paul,
M., Sheta, H., Werner, D. and Ölmann, U., MUFTE-UG: Structure, Applications and
Numerical Methods, Newsletter, International Groundwater Modeling Centre, Col-
orado School of Mine, 2005. Volume XXIII (2), 2005.
Institut für Wasserbau
Universität Stuttgart
Pfaffenwaldring 61
70569 Stuttgart (Vaihingen)
Telefon (0711) 685 - 64717/64749/64752/64679
Telefax (0711) 685 - 67020 o. 64746 o. 64681
E-Mail: [email protected]
http://www.iws.uni-stuttgart.de
2 Marotz, Günter: Beitrag zur Frage der Standfestigkeit von dichten Asphaltbelägen
im Großwasserbau, 1964
7 Sonderheft anläßlich des 65. Geburtstages von Prof. Arthur Röhnisch mit Beiträ-
gen von Benk, Dieter; Breitling, J.; Gurr, Siegfried; Haberhauer, Robert; Hone-
kamp, Hermann; Kuz, Klaus Dieter; Marotz, Günter; Mayer-Vorfelder, Hans-Jörg;
Miller, Rudolf; Plate, Erich J.; Radomski, Helge; Schwarz, Helmut; Vollmer, Ernst;
Wildenhahn, Eberhard; 1967
12 Erbel, Klaus: Ein Beitrag zur Untersuchung der Metamorphose von Mittelgebirgs-
schneedecken unter besonderer Berücksichtigung eines Verfahrens zur Bestim-
mung der thermischen Schneequalität, 1969
15 Schulz, Manfred: Berechnung des räumlichen Erddruckes auf die Wandung kreis-
zylindrischer Körper, 1970
17 Benk, Dieter: Ein Beitrag zum Betrieb und zur Bemessung von Hochwasser-
rückhaltebecken, 1970
Verzeichnis der Mitteilungshefte 3
19 Kuz, Klaus Dieter: Ein Beitrag zur Frage des Einsetzens von Kavitationserschei-
nungen in einer Düsenströmung bei Berücksichtigung der im Wasser gelösten Ga-
se, 1971,
21 Sonderheft zur Eröffnung der neuen Versuchsanstalt des Instituts für Wasserbau
der Universität Stuttgart mit Beiträgen von Brombach, Hansjörg; Dirksen, Wolfram;
Gàl, Attila; Gerlach, Reinhard; Giesecke, Jürgen; Holthoff, Franz-Josef; Kuz, Klaus
Dieter; Marotz, Günter; Minor, Hans-Erwin; Petrikat, Kurt; Röhnisch, Arthur; Rueff,
Helge; Schwarz, Helmut; Vollmer, Ernst; Wildenhahn, Eberhard; 1972
27 Steinlein, Helmut: Die Eliminierung der Schwebstoffe aus Flußwasser zum Zweck
der unterirdischen Wasserspeicherung, gezeigt am Beispiel der Iller, 1972
35 Sonderheft anläßlich des 65. Geburtstages von Prof. Dr.-Ing. Kurt Petrikat mit Bei-
trägen von: Brombach, Hansjörg; Erbel, Klaus; Flinspach, Dieter; Fischer jr., Ri-
chard; Gàl, Attila; Gerlach, Reinhard; Giesecke, Jürgen; Haberhauer, Robert; Haf-
ner Edzard; Hausenblas, Bernhard; Horlacher, Hans-Burkhard; Hutarew, Andreas;
Knoll, Manfred; Krummet, Ralph; Marotz, Günter; Merkle, Theodor; Miller, Chris-
toph; Minor, Hans-Erwin; Neumayer, Hans; Rao, Syamala; Rath, Paul; Rueff, Hel-
ge; Ruppert, Jürgen; Schwarz, Wolfgang; Topal-Gökceli, Mehmet; Vollmer, Ernst;
Wang, Chung-su; Weber, Hans-Georg; 1975
54 Herr, Michael; Herzer, Jörg; Kinzelbach, Wolfgang; Kobus, Helmut; Rinnert, Bernd:
Methoden zur rechnerischen Erfassung und hydraulischen Sanierung von Grund-
wasserkontaminationen, 1983, ISBN 3-921694-54-X
56 Müller, Peter: Transport und selektive Sedimentation von Schwebstoffen bei ge-
stautem Abfluß, 1985, ISBN 3-921694-56-6
69 Huwe, Bernd; van der Ploeg, Rienk R.: Modelle zur Simulation des Stickstoffhaus-
haltes von Standorten mit unterschiedlicher landwirtschaftlicher Nutzung, 1988,
ISBN 3-921694-69-8,
79 Marschall, Paul: Die Ermittlung lokaler Stofffrachten im Grundwasser mit Hilfe von
Einbohrloch-Meßverfahren, 1993, ISBN 3-921694-79-5,
101 Ayros, Edwin: Regionalisierung extremer Abflüsse auf der Grundlage statistischer
Verfahren, 2000, ISBN 3-933761-04-2,
102 Huber, Ralf: Compositional Multiphase Flow and Transport in Heterogeneous Po-
rous Media, 2000, ISBN 3-933761-05-0
108 Schneider, Matthias: Habitat- und Abflussmodellierung für Fließgewässer mit un-
scharfen Berechnungsansätzen, 2001, ISBN 3-933761-11-5
110 Lang, Stefan: Parallele numerische Simulation instätionärer Probleme mit adapti-
ven Methoden auf unstrukturierten Gittern, 2001, ISBN 3-933761-13-1
113 Iqbal, Amin: On the Management and Salinity Control of Drip Irrigation, 2002, ISBN
3-933761-16-6
115 Winkler, Angela: Prozesse des Wärme- und Stofftransports bei der In-situ-
Sanierung mit festen Wärmequellen, 2003, ISBN 3-933761-18-2
120 Paul, Maren: Simulation of Two-Phase Flow in Heterogeneous Poros Media with
Adaptive Methods, 2003, ISBN 3-933761-23-9
121 Ehret, Uwe: Rainfall and Flood Nowcasting in Small Catchments using Weather
Radar, 2003, ISBN 3-933761-24-7
122 Haag, Ingo: Der Sauerstoffhaushalt staugeregelter Flüsse am Beispiel des Ne-
ckars - Analysen, Experimente, Simulationen -, 2003, ISBN 3-933761-25-5
123 Appt, Jochen: Analysis of Basin-Scale Internal Waves in Upper Lake Constance,
2003, ISBN 3-933761-26-3
10 Institut für Wasserbau * Universität Stuttgart * IWS
124 Hrsg.: Schrenk, Volker; Batereau, Katrin; Barczewski, Baldur; Weber, Karolin und
Koschitzky, Hans-Peter: Symposium Ressource Fläche und VEGAS - Statuskol-
loquium 2003, 30. September und 1. Oktober 2003, 2003, ISBN 3-933761-27-1
125 Omar Khalil Ouda: Optimisation of Agricultural Water Use: A Decision Support
System for the Gaza Strip, 2003, ISBN 3-933761-28-0
127 Witt, Oliver: Erosionsstabilität von Gewässersedimenten mit Auswirkung auf den
Stofftransport bei Hochwasser am Beispiel ausgewählter Stauhaltungen des Ober-
rheins, 2004, ISBN 3-933761-30-1
130 Reichenberger, Volker; Helmig, Rainer; Jakobs, Hartmut; Bastian, Peter; Niessner,
Jennifer: Complex Gas-Water Processes in Discrete Fracture-Matrix Systems: Up-
scaling, Mass-Conservative Discretization and Efficient Multilevel Solution, 2004,
ISBN 3-933761-33-6
131 Hrsg.: Barczewski, Baldur; Koschitzky, Hans-Peter; Weber, Karolin; Wege, Ralf:
VEGAS - Statuskolloquium 2004, Tagungsband zur Veranstaltung am 05. Oktober
2004 an der Universität Stuttgart, Campus Stuttgart-Vaihingen, 2004, ISBN 3-
933761-34-4
132 Asie, Kemal Jabir: Finite Volume Models for Multiphase Multicomponent Flow
through Porous Media. 2005, ISBN 3-933761-35-2
133 Jacoub, George: Development of a 2-D Numerical Module for Particulate Con-
taminant Transport in Flood Retention Reservoirs and Impounded Rivers, 2004,
ISBN 3-933761-36-0
134 Nowak, Wolfgang: Geostatistical Methods for the Identification of Flow and Trans-
port Parameters in the Subsurface, 2005, ISBN 3-933761-37-9
135 Süß, Mia: Analysis of the influence of structures and boundaries on flow and
transport processes in fractured porous media, 2005, ISBN 3-933761-38-7
139 Kobayashi, Kenichiro: Optimization Methods for Multiphase Systems in the Sub-
surface - Application to Methane Migration in Coal Mining Areas, 2005,
ISBN 3-933761-42-5
143 Wege, Ralf: Untersuchungs- und Überwachungsmethoden für die Beurteilung na-
türlicher Selbstreinigungsprozesse im Grundwasser, 2005, ISBN 3-933761-46-8
145 Hrsg.: Braun, Jürgen; Koschitzky, Hans-Peter; Müller, Martin: Ressource Unter-
grund: 10 Jahre VEGAS: Forschung und Technologieentwicklung zum Schutz von
Grundwasser und Boden, Tagungsband zur Veranstaltung am 28. und 29. Sep-
tember 2005 an der Universität Stuttgart, Campus Stuttgart-Vaihingen, 2005, ISBN
3-933761-48-4
154 Das, Tapash: The Impact of Spatial Variability of Precipitation on the Predictive
Uncertainty of Hydrological Models, 2006, ISBN 3-933761-57-3
157 Manthey, Sabine: Two-phase flow processes with dynamic effects in porous
media - parameter estimation and simulation, 2007, ISBN 3-933761-61-1
158 Pozos Estrada, Oscar: Investigation on the Effects of Entrained Air in Pipelines,
2007, ISBN 3-933761-62-X
159 Ochs, Steffen Oliver: Steam injection into saturated porous media – process
analysis including experimental and numerical investigations, 2007,
ISBN 3-933761-63-8
160 Marx, Andreas: Einsatz gekoppelter Modelle und Wetterradar zur Abschätzung
von Niederschlagsintensitäten und zur Abflussvorhersage, 2007,
ISBN 3-933761-64-6
162 Kebede Gurmessa, Tesfaye: Numerical Investigation on Flow and Transport Char-
acteristics to Improve Long-Term Simulation of Reservoir Sedimentation, 2007,
ISBN 3-933761-66-2
166 Freeman, Beau: Modernization Criteria Assessment for Water Resources Plan-
ning; Klamath Irrigation Project, U.S., 2008, ISBN 3-933761-70-0
Verzeichnis der Mitteilungshefte 13
168 Yang, Wei: Discrete-Continuous Downscaling Model for Generating Daily Precipi-
tation Time Series, 2008, ISBN 3-933761-72-7
173 Wagner, Sven: Water Balance in a Poorly Gauged Basin in West Africa Using At-
mospheric Modelling and Remote Sensing Information, 2008,
ISBN 978-3-933761-77-4
174 Hrsg.: Braun, Jürgen; Koschitzky, Hans-Peter; Stuhrmann, Matthias; Schrenk, Vol-
ker: VEGAS-Kolloquium 2008 Ressource Fläche III, Tagungsband zur Veranstal-
tung am 01. Oktober 2008 an der Universität Stuttgart, Campus Stuttgart-
Vaihingen, 2008, ISBN 978-3-933761-78-1
175 Patil, Sachin: Regionalization of an Event Based Nash Cascade Model for Flood
Predictions in Ungauged Basins, 2008, ISBN 978-3-933761-79-8
179 Laux, Patrick: Statistical Modeling of Precipitation for Agricultural Planning in the
Volta Basin of West Africa, 2009, ISBN 978-3-933761-83-5
180 Ehsan, Saqib: Evaluation of Life Safety Risks Related to Severe Flooding, 2009,
ISBN 978-3-933761-84-2
Die Mitteilungshefte ab der Nr. 134 (Jg. 2005) stehen als pdf-Datei über die Homepage
des Instituts: www.iws.uni-stuttgart.de zur Verfügung.