جامعة تونس المنار [301438]
جامعة تونس المنار
Université de Tunis El Manar
THESE
présentée pour obtenir le titre de
DOCTEUR DE L’ECOLE NATIONALE D’INGENIEURS DE TUNIS
École doctorale : Sciences et Techniques de l’Ingénieur
Spécialité : Génie Hydraulique
Par
EZZEDDINE LAABIDI
REACTIVE TRANSPORT MODELING IN SATURATED POROUS MEDIA: [anonimizat], water-[anonimizat], pore and particle size distributions. Porosity and permeability can be changed due to dissolution and precipitation of minerals. [anonimizat] ([anonimizat]-freshwater mixing zone). [anonimizat] a [anonimizat], investigative tools and management of such systems. The majority of the studies quantifies these hydrogeological systems using traditional density dependent flow and transport models and do not consider the effect of the chemical reactions. Interdependence of density dependent flow and chemical reactions and their effects on the porosity and permeability are the most important keys towards a reliable modeling of these complex systems. The dissolution or the precipitation of such rocks can easily induce a development of the porosity and permeability as a result of the mixing processes in the coastal carbonate aquifers (Laabidi and Bouhlila 2015; Romanov and Dreybrodt 2006). The modeling of such problem requires a set of highly nonlinearly coupled equations. GEODENS code (Bouhlila 1999b; Bouhlila and Laabidi 2008a; Laabidi and Bouhlila 2015), [anonimizat] a finite element procedure. [anonimizat]. Its main purpose is to represent the physicochemical processes in the subsurface system. [anonimizat]. GEODENS (Bouhlila 1999b; Bouhlila and Laabidi 2008a; Laabidi and Bouhlila 2015) comprises two modules, (i) the density dependent flow and multispecies transport module and (ii) the geochemical module. The mathematical formulation of the first module leads to a nonlinear and strongly coupled equations and in order to solve such equations, a finite element method has been developed with a consistent numerical scheme of gravity terms to calculate Darcy velocities (C. I. Voss 1984b; C. I. Voss and Souza 1987). [anonimizat] (Pitzer et al. 1984). This geochemical module allows the calculation of ions and solvent activities as well as the density of the solution. [anonimizat]- precipitation processes are controlled by diffusion and are represented by a first order kinetic law. [anonimizat]. The first one focuses on the effect of precipitation reaction on the flow and transport problem in a 2D flow cell (Laabidi and Bouhlila 2016a). [anonimizat] and precipitation continue along the centerline cell (the mixing line). The effect of calcite deposition on the flow and transport is discussed in term of two parameters: porosity decrease and the total aqueous calcium concentrations. The second simulation treated the effect of calcite dissolution on the seawater intrusion in coastal carbonate homogeneous aquifer (Laabidi and Bouhlila 2016b). This problem is also modeled using a semi-analytical method (Laabidi and Bouhlila 2015). Through the results of this simulation, we have shown that the change of porosity is created in narrow regions at the freshwater side of the mixing zone. The behavior of the solution is discussed in terms of two parameters: velocity field which has a direct effect on the seawater flux, and the penetration length. An increase of the penetration length of the seawater is identified during calcite dissolution. Finally, it is hoped that the code will be a useful tool for understanding hydrological and geochemical processes in arid regions.
Keywords: Finite element method, variable density flow, calcite precipitation reaction; clogging effect; kinetic precipitation reaction; porosity change, seawater intrusion, calcite dissolution reaction; reactive Henry problem; multispecies reactive transport; kinetic dissolution reaction; penetration length; coastal carbonate aquifer.
Résumé
Dans plusieurs situations hydrogéochimiques, les espèces carbonatées interagissent avec la matrice solide selon des processus de dissolution et de précipitation : (1) Lors de l’évaporation des saumures évoluant dans la voie carbonatée, les premiers sels déposés sont calciques ; calcite et dolomite. Ils constituent souvent les premières couches des dépôts évaporitiques des Sebkhas, Chotts et autres Playas. (2) Dans un contexte de réchauffement climatique, la séquestration du CO2 dans les formations géologiques profondes est une solution de plus en plus évoquée. Il est clair que l’injection de ce gaz dans des roches calcaires va induire des processus de dissolution ce qui devra à plus ou moins long terme compromettre le confinement de ces pièges et la fiabilité de cette solution. (3) Dans les aquifères côtiers carbonatés, l’intrusion marine peut être aggravée du fait de la dissolution des roches calcaires dans le mélange eau douce – eau de mer. La quantification de ces différents processus est nécessaire en vue de prédire leurs impacts sur les divers environnements en question. La modélisation hydrogéochimique couplée de ces phénomènes dans les milieux poreux constitue un outil d’analyse et de prédiction de ces processus. Dans le cas de mélange des eaux salées et des eaux douces dans les aquifères carbonatés, objet de ce travail, il s’agit de préciser les variations de la porosité suite à la dissolution ou la précipitation de la calcite, de quantifier la variation conséquente des perméabilités et ce qui s’ensuit comme modifications du champ de l’écoulement et des concentrations en solutés dans les franges côtières de ces nappes. Le plus souvent, ces processus de dissolution ou de précipitation, conjugués avec la surexploitation des ressources, constituent un risque réel pour la qualité des eaux dans ces systèmes. Nous présentons dans ce travail GEODENS un modèle hydrogéochimique couplé qui simule ces processus (Bouhlila 1999b; Bouhlila and Laabidi 2008a; Laabidi and Bouhlila 2015). L’écoulement densitaire et le transport multi espèces couplés y sont modélisés par une méthode aux éléments finis avec une procédure de pondération mixte des termes de gravité dans l’équation de Darcy, capable de rendre compte des courant de convection densitaire (C. I. Voss 1984b; C. I. Voss and Souza 1987). L’état de saturation de la solution aqueuse vis-à-vis des sels évaporitiques est testé en chaque point et à chaque instant dans le domaine d’étude et les progressions des réactions de dissolution ou de précipitation sont alors écrites selon une description cinétique. Les volumes libérés ou occupés par les sels dissous ou précipités sont algébriquement quantifiés dans le volume de pores et les perméabilités sont modifiées en conséquence selon des expressions basées sur des relations empiriques. Du fait des forces ioniques élevées dans le domaine des sels et des saumures qui nous concerne, nous avons retenu le modèle de Pitzer (Pitzer et al. 1984) pour le calcul des activités chimiques des espèces aqueuses et de l’eau. Dans ce travail, le code GEODENS vérifié et validé est utilisé pour la modélisation de l’impact des réactions de dissolution et de précipitation de la calcite dans les zones de mélanges sur les écoulements et le transport dans les aquifères calcaires. La première simulation porte sur l'effet de la réaction de précipitation sur l'écoulement et le transport dans une cellule d'écoulement 2D (Laabidi and Bouhlila 2016a). Le principal résultat montre que la précipitation de la calcite peut provoquer le colmatage de l’aquifère au niveau de la ligne de mélange. L'effet du dépôt de la calcite sur l’écoulement et le transport est discuté en terme de deux paramètres: la diminution de la porosité et la concentration totale de calcium. La deuxième simulation traite l'effet de dissolution de la calcite sur l'intrusion d'eau de mer dans des aquifères de type calcaire (Laabidi and Bouhlila 2016b). Ce problème est traité aussi par une méthode semi-analytique (Laabidi and Bouhlila 2015). Le principal résultat de cette modélisation montre que ce problème est, contrairement au problème de Henry classique, instationnaire et que le degré d’intrusion des eaux de mer dans les nappes côtières est en progression perpétuelle. Enfin, nous souhaitons que le code GEODENS soit un outil utile pour la compréhension des processus hydrologiques et géochimiques dans les régions arides.
Mots-clés: méthode des éléments finis, écoulement densitaire, la réaction de précipitation de calcite; effet de colmatage; cinétique des réactions de précipitation; changement de porosité, intrusion marine, réaction de dissolution de calcite; problème d’Henry réactif; transport multiespèces réactif; cinétique des réactions de dissolution; longueur de pénétration; aquifère côtière carbonaté.
CONTENT
CHAPTER 1. INTRODUCTION 5
I.II Mixing inducing calcite dissolution 6
I.II Mixing inducing calcite precipitation 7
I.III Aims and structure of thesis 9
CHAPTER 2. LITTERATURE OVERVIEW AND PREVIOUS MODEL DEVELOPMENT 12
II.I Carbonate rock 13
II.I.2. Carbonate diagenesis 13
II.I.2.1.Carbonate ramp and carbonate reef shelf 13
II.I.2.2. Classification of Carbonate Porosity 15
II.II Carbonate Hydrogeology 16
II.III. Coastal aquifer and density dependent flow model 20
II.III.1 The most known situation of seawater intrusion 21
II.III.1.1 Horizontal intrusion 21
II.III.1.2 Downward leakage of saltwater from surface water 21
II.III.1.3 Upconing below pumping wells 22
II.III.2 Survey of Numerical model and benchmarking density dependent flow models: 23
II.IV Geochemical Settings 27
II.IV.1 Kinetic reactions in sedimentary basin 29
II.IV.2 Mixing of carbonate waters 31
II.V Reactive transport modeling in variable density flow 34
II.VI Reactive transport modeling in carbonate coastal aquifer 40
CHAPTER 3. GOVERNING EQUATIONS 45
III.I. Density-dependent flow and transport model 45
III.II. Geochemical Model 48
III.III. Hydrogeochemical model: Coupled model 51
III.III.1. Governing equations 51
III.III.2. Existing permeability–porosity relationships and implementation into reactive transport models 52
CHAPTER 4. NUMERICAL INTEGRATION 54
IV.I. Integration of coupled equations 54
IV.II. Calculation strategy and resolution 59
CHAPTER 5. MODEL VERIFICATION 63
V.I. Benchmarking Density dependent flow model 63
V.I.1. Rayleigh problem 63
V.I.2. Henry problem 69
V.II. Benchmarking geochemical model 71
V.III. Benchmarking the coupled model 73
V.III.1. Dissolution experiment 75
V.III.2. Precipitation experiment 77
CHAPTER 6. MODEL APPLICATION AND ALTERNATIVE APPROACH 81
VI.I Model Application 81
VI.I.1. Precipitation reaction simulation 82
VI.I.2. Dissolution reaction simulation: Reactive Henry problem 90
VI.II Alternative modeling approach 104
VI.II.1. Mixing of carbonate waters 106
VI.II.2.Alternative modeling approach: mathematical formulation and protocol of the numerical procedure 109
VI.II.2.1. Governing equation 109
VI.II.2.2. Implementation of porosity-permeability relationship 112
VI.II.2.3 Protocol used for the calculation 113
VI.II.3. Testing the approach 115
VI.II.4. Simulation and discussion : Henry Reactive problem: 117
CHAPTER 7. CONCLUSIONS AND PERSPECTIVES 123
VII.I Mixing inducing calcite dissolution 124
VII.I.1. Using the coupled hydrogeochemical model 124
VII.I.2. Using anlutical solution 125
VII.II Mixing inducing calcite precipitation 125
VII.III Perspectives 126
APPENDIX 128
REFERENCES 136
List of figures
Figure II.4. . 14
Figure II.5. . 19
Figure II.6 22
Figure II.7 31
Figure II.8 33
Figure IV.1 55
Figure IV.2 55
Figure IV.3 62
Figure V.1 65
Figure V.2 68
Figure V.3 69
Figure V.4 70
Figure V.5 71
Figure V.6 72
Figure V.7 74
Figure V.8. 75
Figure V.9 76
Figure V.10 78
Figure V.11 79
Figure.VI.1. . 84
Figure.VI.2.. 85
Figure.VI.3.. 86
Figure.VI.4. 88
Figure.VI. 5. 89
Figure.VI.6. 90
Figure.VI.7 93
Figure.VI.8 94
Figure.VI.9 95
Figure.VI.10 96
Figure.VI.11 97
Figure.VI.12 98
Figure.VI.13 100
Figure.VI.14 . 103
Figure.VI.15 104
Figure.VI.16. 107
Figure.VI.17. 108
Figure.VI.18. 114
Figure.VI.19. 115
Figure.VI.20. 116
Figure.VI.21. 117
Figure.VI.22. 118
Figure. VI.23. 119
Figure. VI.24. 120
Figure.VI.25 122
List of tables
Table II.1 17
Table II.2 26
Table II.3 39
Table V.1 65
Table V.2 67
Table V.3 70
Table V.4 72
Table V.5 74
Table V.6 77
Table V.7 80
Table VI.1. 84
Table VI.2. 91
Table VI.3 92
Table VI.4 92
Table VI.5. 101
Table VI.6 121
Table VI.7 121
CHAPTER 1. INTRODUCTION
Environmental engineering faces serious challenges with the management of groundwater contamination, waste disposal, radioactive disposal and greenhouse gas emissions due to anthropogenic activities and their consequences (Pauline et al. 2011; Rechard 1999). In the last few decades interest is increasingly being focused upon predicting how hydrogeochemical systems will evolve over long periods of time. Modeling and computer simulation is a valuable tool that can be used for the understanding of such processes both to interpret laboratory experiments and field data as well as to make predictions and scenarios of long term behavior. One of the most recognized reactive transport problems is the dissolution and precipitation reaction during mixing of seawater and freshwater in coastal carbonate aquifer. An important aspect of reactive transport modeling is the appropriate evaluation of porosity and permeability development during dissolution reaction and their impact on flow and transport. In the following two situations involving these processes are cited:
From the beginning of the twentieth century, underground storage in salty layers is used. Currently, salty formation in deep geological layers is used for the burial of industrial long-life waste, including radioactive waste and CO2 in order to counteract the greenhouse effect. Suitable and safe storage sites are searched and costly investigated in many places around the word (Ghanbari et al. 2006; Hellevang et al. 2005; Kang et al. 2010). The injection of CO2 in the subsurface leads to a relatively rapid decrease of the water solution pH. This decrease will then be buffered by mineral reactions, such as calcite dissolution (Lagneau et al. 2005). Calcite dissolution can induce an increase in the porosity and can enhance further leakage and provide more pathways for the CO2 (Caldeira and Rau 2000).
Soil and groundwater salinization due to intense evapotranspiration in arid regions and/or saline water intrusion from the sea or from Chotts or Sebkhas (playas) is a limiting factor for development and food security in these areas (M. Tirado et al. 2010). However, the relevance of Sebkhas and Chotts in terms of mineral resources is important (Bouhlila 1999b). The Tunisian territory, due to its hydroclimatic conditions, includes a large number of such depressions which contain many salts. The exploitation of these deposits cannot be done without a good knowledge of these environments, as well as of the mechanisms of formation and renewal of salts and brines (Y Yechieli and Wood 2002; Nasri et al. 2015).
The common feature of all these issues, and obviously many others related to salts and brines in porous media is the coupling between geochemical processes, mainly dissolution and precipitation reaction, and physical ones of density-dependent flow and multi-species reactive transport. In order to model these hydrogeological systems, a general model that can handle density dependent flow and multispecies reactive transport modeling is required. The prediction of mineral interactions and water interactions using mathematical models is integral to the study of the quality and the geochemistry of natural water.
I.II Mixing inducing calcite dissolution
Freshwater stored in coastal aquifer is extensively extracted through pumping wells due to high water demand in coastal area (touristic, industrial and public use). In order to enhance freshwater security and to avoid contamination of water reserves, seawater intrusion becomes a topic of great interest for hydrogeologist. During the last few decades, hydrogeologist provided a deeper understanding of the prediction, processes, investigative tools and management of such systems. The majority of the studies quantifies these hydrogeological systems using traditional density dependent flow and transport models and do not consider the effect of the chemical reactions. Interdependence of density dependent flow and chemical reactions and their effects on the porosity and permeability are the most important keys towards a reliable modeling of these complex systems. Seawater intrusion can increased by the solid matrix dissolution processes in the saltwater–freshwater mixing zone in the case of coastal carbonate aquifer. Many investigations have been published in the last few years regarding the hydrogeochemical modeling, especially the issues related to calcite dissolution reaction during mixing of two different waters. The most known studies are as follow:
(Sanford and Konikow 1989b) provided a quantitative estimate of the amount of calcite dissolution expected in mixing zone under typical hydrodynamic and geochemical conditions. To assess the effects of calcite dissolution on the rock matrix, the authors proposed a procedure to quantify the effect on the rock proprieties (porosity and permeability). They used the geochemical code PHREEQC (Parkhurst and Apello 1999b) coupled with a finite difference model for variable density groundwater flow and solute transport.
(Rezaei et al. 2005) extended the study of (Sanford and Konikow 1989b) to simulate the porosity development in coastal aquifer due to calcite dissolution. The authors used RETRASO code, a finite element model (Maarten W Saaltink et al. 2004) to solve simultaneously the aqueous complexation reaction, the dissolution-precipitation reaction and the flow and transport problem using a global implicit method. They found through the simulations two locations of maximum dissolution as obtained by (Sanford and Konikow 1989b), one near the aquifer bottom and another near the discharge boundary.
(Romanov and Dreybrodt 2006) used the analytical solution of (Phillips 1991) to set up an alternative modeling approach to determinate the porosity change in the saltwater-freshwater mixing zone of coastal carbonate aquifers. As a first step, they solve the advection transport equation for salinity s using the finite difference code SEAWAT (Guo and Langevin 2002b). The second step consists in calculating the dissolution capacity in the mixing zone by estimating the calcium equilibrium concentration that is obtained as a function of salinity using PHREEQC (Parkhurst and Apello 1999b). They applied this approach for a stationary salinity distribution.
(Laabidi and Bouhlila 2015) used the analytical solution of (Phillips 1991) in non-stationary step by step regime. At each time step the quantity of the dissolved calcite is quantified, the change of porosity is calculated and the permeability is updated. The reaction rate which is the second derivate of the calcium equilibrium concentration in the equation is calculated by using PHREEQC code (Parkhurst and Apello 1999b). The authors used this result in a 2D finite element code GEODENS (Bouhlila 1999b; Bouhlila and Laabidi 2008a; Laabidi and Bouhlila 2015) to calculate the initial change of the porosity after calculating the salinity gradient. For the next time step the same protocol is used but by using the updated porosity and permeability distributions.
I.II Mixing inducing calcite precipitation
For a given mineral, if the ion product is greater than the solubility product, then the solution is supersaturated with respect to this mineral and precipitation can occur (W. Stumm and Morgan 1970). Supersaturation may be occurred by three main mechanisms: (i) changes of fluid pressure or/and temperature during flow through porous media, (ii) dissolution of a given mineral that results in a supersaturation with respect to another one (e.g. dissolution of calcite can induce precipitation of gypsm (Singurindy and Berkowitz 2003a)), and (iii) mixing induced precipitation (Chen et al. 2013; Hui and Peiyue 2011, 2012; Yoon et al. 2012). Referring to (L. N. Plummer 1975a) the two most important factors generating the precipitation or the dissolution of such solids which are the amount of the seawater in the mixture and the pCO2. If the fraction of the seawater in the mixture zone is high and pCO2 is low then, the precipitation of the calcite in the groundwater can be generated. These factors could explain the absence of dissolution in the coastal aquifer of Mallorca, Spain (Price and Herman 1991a).
Mineral deposition in geological formation and decrease of porosity due to the infilling of fractured and porous media during precipitation reactions are of great interest in many environmental and industrial applications such as effluent disposal (Eriksson and Destouni 1997) and oil recovery extraction (Rege and Fogler 1989). The lack of data on precipitation kinetics in subsurface is the main cause of the miss understanding of calcite filling in porous media. Modeling of expected time and amount of calcite needed to produce mineral filling has been a largely unresolved problem. In this context a few works have been carried out toward understanding reactive transport modeling associated with porosity change during precipitation reaction and in the following we quote some of these works:
(Lee and Morse 1999) developed a new experimental protocol in order to investigate the precipitation patterns of calcite along a synthetic vein channel, to study the filling processes of calcite, and to estimate the fluid volume and the time needed to produce natural calcite veins from experimental and simulation results. During mixing, precipitation can occur where region of high porosity can be filled while other can remain unchangeable (Emmanuel and Berkowitz 2005).
(Hilgers and Urai 2002; Lee and Morse 1999) found that when supersaturated fluids flowed through artificial fractures, most of the mineral deposition occurred within several centimeters of the inlet. Furthermore, (Hilgers and Urai 2002) reported that by the end of the experiment, the inlet had become completely clogged and deposition started ‘‘downstream’’.
(Emmanuel and Berkowitz 2005) conducted a series of simulations in order to determinate the porosity evolution during precipitation of calcite in a 2D modeled domain using a finite element scheme. The authors checked the dynamics and patterns of porosity change during mixing induced calcite precipitation in a geological realistic system (homogeneous and heterogeneous system containing high porosity region). The authors also examined the specific surface area using two expressions function of porosity. Their model represents comprehensive tools for calcite precipitation in porous and fractured media (mixing induced precipitation can account for system in which some high porosity region are clogged while other remains unchangeable).
(Singurindy and Berkowitz 2003a) performed a laboratory experiment to understand the evolution of hydraulic conductivity and flow patterns caused by dissolution of dolomite and precipitation of calcium carbonate in porous dolomite rock by using a column flow experiments. They conclude that long term experiments show that the evolution of hydraulic conductivity in carbonate rock, due to competition among flow, precipitation and dissolution processes, can display a wide variety of possible behaviors.
(Tartakovsky et al. 2008) conducted a series of laboratory experiments in order to investigate pore scale and continuum scale mixing inducing precipitation. Two solutions containing Na2CO3 and CaCl2 were injected in parallel in a 2D flow cell filled with quartz sand. Several simulations were carried out by the authors for a wide range of Peclet and Damköhler numbers. The authors found that calcite precipitate near the inlet boundary in the mixing zone between the two solutions. The authors also found that during precipitation, clogging phenomenon can be generated. This phenomenon separates the two solutions, then they may become undersaturated with respect to calcite and dissolution can occur.
(Katz et al. 2011) conducted a laboratory experiment in order to provide a quantitative analysis of calcite evolution patterns during mixing of two solutions of Na2CO3 and CaCl2 equilibrated with calcite. A 2D flow cell with dimension of 25cmx10cmx0.8cm in the x-, y- and z-direction respectively is used. The authors found that calcium carbonate precipitated along the mixing zone between the two solutions which flowed in parallel. The precipitation of calcium carbonate causes a decrease of porosity as pores in the mixing zone are filled by the precipitated mineral. In order to evaluate the effect of porosity change on the flow and transport problem, they measured the Ca2+ concentration in some selected sampling ports and used the RETRASO code (M. W. Saaltink et al. 2004) to numerically reproduce the results found in laboratory experiment. They found that the model appear to capture the main features of the processes of the evolution of calcium concentration in flow domain.
I.III Aims and structure of thesis
The objective of the thesis is to present a general numerical tool for a multicomponent reactive transport. The model was developed by (Bouhlila 1999b) and has the capability to simulate physical and geochemical processes in porous media in one and two dimensions. In order to simulate calcite dissolution and precipitation reaction during mixing of seawater and freshwater, GEODENS is extended to carbnates species (Bouhlila and Laabidi 2008a). To achieve this goal it is necessary to develop a protocol to couple geochemical reactions (mainly salt dissolution and precipitation reaction), flow and transport problem. In this regard, the dissolved or precipitated salt need to be quantified, the change of porosity and permeability need to be calculated and the effect of this later on the reactive transport is evaluated.
The present thesis is divided into five principle chapters presented in submitted and accepted article. In the second chapter, we present a literature overview of the main futures of carbonate coastal aquifer. The objectives of the second chapter are (1) to describe the geological and hydrogeological settings of carbonate aquifer, (2) to describe the density dependent flow in coastal aquifer and the geochemical setting of mixed carbonate waters, and (3) to survey a literature overview of the previous reactive modeling study in density dependent flow and carbonate coastal aquifer.
In the chapter 3, we present the mathematical framework of the GEODENS code used in this work. The fundamental mathematical relationships that are used to describe the relevant processes are introduced. GEODENS code used in this work comprises two modules, (i) the density dependent flow and multispecies transport module and (ii) the geochemical module. The mathematical formulation of the first module leads to a nonlinear and strongly coupled equations. The second module focuses on salts and brine geochemistry using the Pitzer model. This geochemical module allows the calculation of ions and solvent activities as well as the density of the solution. Reactions of salts including dissolution- precipitation processes are controlled by diffusion and are represented by a first order kinetic law.
The governing equations presented in chapter 3 are discretized in chapter 4. This chapter introduces the numerical techniques used for the resolution of the governing equations for variably-saturated flow and reactive transport. The governing equations are approximated by discretized equations using a finite element method. By the application of Galerkine method and by using integration by parts we find the matrix form of the reactive transport problem and the discontinuity velocity approach is used to evaluate the Darcy velocities. Finally, the time discretization scheme used in this work is implicit.
Chapter 5 presents the verification of GEODENS based on five test problems. Regarding the density dependent flow model, the first benchmark used is the Henry problem, (H. R. Henry 1964a) solved the governing equations of the density-dependent flow in a non-dimensional form and presented a semi-analytical solution for the steady state distribution of salt concentration and stream functions. The second benchmark used is the convection cell problem which describes the free-convection induced by buoyancy effects. The geochemical component of the GEODENS code was successfully verified using results on the effect of the saltwater fraction on the chemical processes (dissolution-precipitation reactions) in a carbonate environment during the mixing between saltwater and freshwater, presented by (Rezaei et al. 2005). And finally, the coupled model was validated by using the laboratory presented by (Singurindy et al. 2004a), a simulations of the porosity evolution in a 2D flow cell due to calcite precipitation-dissolution reactions was carried out and compared to the laboratory experiments results.
Chapter 6 presents the applicability of GEODENS for the simulation of reactive transport models. The first part of this chapter focuses on the capability of GEODENS code (Bouhlila 1999b; Bouhlila and Laabidi 2008a; Laabidi and Bouhlila 2015) in calcite dissolution and precipitation reactions and their effect on the flow and transport problem. The first application treated the clogging effect on the flow and transport problem during calcite precipitation reaction and the second application treated the effect of the calcite dissolution on the penetration length of the seawater intrusion. The second part of this chapter deals with the Question: the numerical coupled model is it the unique solution of reactive transport modeling problem? In this regard, we present in the second part of this chapter an alternative modeling approach to solve reactive transport modeling during calcite dissolution and precipitation reactions based on the analytical solution presented by (Phillips 1991).
The seventh and last chapter summarizes and discusses the principal results. In this conclusion we are also proposing a perspective of eventual future fields of research to carry on as a continuation of the present work.
CHAPTER 2. LITTERATURE OVERVIEW AND PREVIOUS MODEL DEVELOPMENT
Carbonate are sedimentary rocks containing more than 50% carbonate minerals which are mainly calcite, CaCO3, and dolomite, CaMg(CO3)2. The carbonate minerals could be a result of chemical precipitation, organic processes or may occur as detrital material. Some of the dolomites may be a result of diagenesis known as dolomitization which involves replacement of calcite by dolomite to a varying extent. The term limestone is used for those rocks which contain more than 90% of carbonates. If the rock contains more than 50% but less than 90% carbonates, it is termed as arenaceous limestone or an argillaceous limestone, depending upon the relative amounts of quartz and clay minerals. Chalk is a type of limestone which is soft and white in colour and is rich in shell fragments. The carbonate rocks occupy about 10% of the Earth’ surface and supply nearly a quarter of the world’s population with water. Locally, their cumulative thickness may be 10 000 m (Singhal and Gupta 2010). In the geological past, most limestones and dolomites were deposited during the intermediate phases of Caledonian, Hercynian and Alpine tectonic cycles and their ages are mainly Palaeozoic and Mesozoic. One of the largest carbonate aquifers in the world is the Floridan aquifer system of Palaeocene to Miocene age in the southeast United States consisting of gently dipping thick sequences of carbonate sediments separated by less permeable clastic sediments. Carbonate rocks have attracted greater attention of hydrogeologists in the past few decades, both from the point of view of water supply as well as potential sites for waste disposal. Carbonate rocks are also of great economic value having large deposits of oil, gas and minerals. Carbonate reservoirs often represent a source of high quality water and containing more than 60% of the word’s oil reserves and 40% of the word’s gas reserves (Nurmi and Standen 1997). Carbonate reservoirs representing suitable hosts in underground storage and the high reactivity of the rock causing significant changes in their structure with the creation of barriers permeability or to the development of dissolution zones during mixing of freshwater and seawater (cavities, karsts) that may affect the water resources in the case of an aquifer (Purser 1980).
The objectives of this chapter are (1) to describe the geological and hydrogeological settings of carbonate aquifer, (2) to describe the density dependent flow in coastal aquifer and the geochemical setting of mixed carbonate waters, and (3) to survey a literature overview of the previous reactive modeling study in density dependent flow and carbonate coastal aquifer.
II.I Carbonate rock
II.I.2. Carbonate diagenesis
Carbonate rocks have a strong diagenetic potential. A variations of temperature and pressure or chemical composition of the fluid generate rapid changes in mineralogy and porosity of the rock. Also it should be mentioned that carbonates are susceptible to rapid and extensive diagenetic change. Carbonate minerals are susceptible to rapid dissolution, cementation, recrystallization, and replacement at ambient conditions in a variety of diagenetic environments. Porosity and permeability in carbonate reservoirs depend on a broad array of rock properties, on diagenetic episodes that may continue from just after deposition through deep burial, and on fracture patterns related more to the geometry of stress fields than to rock type (Lucia 1999).
II.I.2.1.Carbonate ramp and carbonate reef shelf
Carbonates, like sandstones, occur in identifiable sequences which reflect changing marine conditions and environments (Fig. II.4). Carbonates can be deposited in a wide range of marine environments. They typically occur in sequences which can be characterized as ramp (a) or reef shelf (b) settings. Low-energy environments, such as the back reef shoals, which are protected from wave and current action, are characterized by higher concentrations of lime mud while clean rocks with high original permeabilities are found in high-energy zones at the shoreline or around the main reef wall.
The term carbonate ramp was adopted to describe a gently sloping depositional surface which passes gradually without slope break from a shallow, high-energy environment to a deeper, low-energy environment. A ramp is attached to a shoreline at one end and to cor- relative basinal beds at the other end. Fig. II.4 exhibits the basic differences between immed and unrimmed shelves (platforms) and ramps with respect to the position of the optimum carbonate production (carbonate factory) and sediment transport.
Carbonate shelf' characterizes a carbonate depositional system which develops a constructional relief above the sea floor and is bordered on the shoreward side by marginal-marine or continental sediments, and offshore by slope and basinal sediments. The transition from shallow water to basin occurs over a relatively short distance and is marked by a distinct break in the slope. In this definition the term means the same as 'attached platform' a term, sometimes used for platform settings starting at the coast and continuing offshore.
Figure II.4. Carbonate ramp and reef shelf (Nurmi and Standen 1997).
II.I.2.2. Classification of Carbonate Porosity
Carbonates have much more complex pore system as compared to the siliciclastics. This complexity is due to biological origin of carbonates and chemical reactivity during and at later stage of deposition. Biological process results in porosity within grains, fossil-related shelter porosity and growth framework porosity in reefs while chemical reactivity due to diagenetic processes, i.e, solution and dolomitization, result in secondary porosity. There are certain other intermediate processes by which porosity is developed, reduced, destructed and modified throughout the burial history.
II.I.2.2.1 Primary Porosity
In depositional reservoirs porosity is formed only by depositional processes and this porosity is classified as primary porosity. Since carbonate rocks are detrital, biogenic or chemical origin. Primary porosity is any porosity present in a sediment or rock at the termination of depositional processes. Primary porosity is formed in two basic stages, the predepositional stage and the depositional stage. The predepositional stage begins when individual sedimentary particles form and includes intragranular porosity such as is seen in forams, pellets, ooids, and other nonskeletal grains. This type of porosity can be very important in certain sediments. The depositional stage is the time involved in final deposition, at the site of final burial of sediment or a growing organic framework. Porosity formed during this stage is called depositional porosity and is important relative to the total volume of carbonate porosity observed in carbonate rocks and sediments (C. H. Moore 1989).
II.I.2.2.2 Secondary Porosity
Secondary porosity is generated and modified due to processes such as dissolution and dolomitization. Fracturing increases the permeability rather than increasing the porosity (Lucia 1999). diagenesis refers to all the post depositional but before metamorphism changes to the sedimentary rocks. All changes in grain size, shape, volume, porosity, chemical composition and change in sedimentary structures happens after the deposition. diagenesis may be of mechanical, biological, chemical origin or combination of these ones. Mechanical diagenesis through compaction might result in volume reduction, interstitial water expulsion and modification in grain packing thus changes and modifies the rock properties after deposition(C. H. Moore 1989) .
II.I.2.2.3 Karst
Due to dissolution of carbonate rocks in carbonic rich water form karsts. Numerous factors play an important role in the development of karst including lithology, hydraulic conductivity, chemically enriched active ground water, palaeoclimates, tectonic activity and eustatic changes etc. Carbon cycle plays an important role in the development of karsts and associated features. Change in hydraulic gradient due to sea level change will favour solution activity and consequently formation of karsts and karst related features. The direction of ground water movement, permeability, lithology and overburden rock thickness bear an important relation with karstification. In vadose zone water moment through fractures will form vertically oriented pipes and sinks while its horizontal movement in phreatic zone will result in horizontally oriented cavities and cave system (Loucks 1999). Exposure of carbonate rocks and unconfined carbonate aquifer, due to recharge and discharge, will tend to enhance the karstification. Dissolution is a dominant active process for karstification in vadose and phreatic zones (C. H. Moore 1989).
II.II Carbonate Hydrogeology
Carbonate aquifer is classified as one of the most important water supply sources and have been the subject of geochemical investigations since the 1970s. As a result, geochemical processes in carbonate aquifer systems and their hydrogeological implications are now well documented for many different sites around the word (Capaccioni et al. 2001; Grasby and Betcher 2002; Hanshaw and Back 1979; Hess and White 1993; Jacobson and Wasserburg 2005; Kohfahl et al. 2008; McIntosh and Walter 2006; L. Niel Plummer et al. 1990; Stetzenbach et al. 1999; Wang et al. 2006; Wicks and Herman 1996).
The application of any groundwater evaluation techniques is satisfactory as long as the hydrogeological settings of a test site satisfy some assumptions. A general assumption that must be checked is the Darcian flow (Hickey 1984). Carbonate aquifers have long represent a complex system for hydogeologist because of the complexity of groundwater flow where it is not trivial to verify the Darcian flow hypothesis (W. B. White 1969). In order to recognize the behavior of carbonate aquifer it is important to understand the groundwater flow in such system. (W. B. White 1969) had developed a useful classification of carbonate aquifer based on geological and hydrogeological settings. The classification is established between (1) diffuse flow, (2) free flow, and (3) confined flow as shown in table II.1.
Table II.1 Types of carbonate aquifer systems in regions of low to moderate relief (W. B. White 1969).
Diffuse flow aquifer: where groundwater flow has been retarded by lithological factors. Aquifer showing fine textured permeability (Legrand and Stringfield 1971). The fine textured permeability may be a result of uniformly spaced and poorly developed fracture systems. As shown in Fig. II.5 (a) cavities are limited in number and size. True caves are rare, small and poorly integrated. A known example is founded in the Silurian dolomite aquifer of the DuPage County-Chicago Region of Illinois (Sasman 1974).
Free flow aquifer: As shown in Fig. II.5 (b), free flow is localized into integrated system of conduits where permeability is developed because of the non-uniform distribution of fracture system that intersects the carbonate body. Large flow takes place in conduits and channels where flow velocity is in turbulent regime and may be in the order of feet per second (Shuster and White 1971). Free flow as defined by (W. B. White 1969) has been subdivided into several categories (see table II.1).
Figure II.5. Schematic diagram of diffuse flow carbonate aquifer (a), free flow aquifer system (b), artesian carbonate aquifer system (c) and sandwich carbonate aquifer system (d) .
Confined flow aquifer: depending on the geological settings of the aquifer, confined flow is subdivided into two subvarieties:
Artesian aquifers (Fig. II.5 (c)): are created by capping beds which are tiled in such a way as to force the groundwater flow at depth under hydrostatic head condition (Legrand and Stringfield 1971).
Sandwich aquifers (Fig. II.5 (d)): are the extreme limit of carbonate unit which is both capped and perched, and is developed by a thin beds of soluble rocks compared to the total thickness of beds above base level. Flow is retarded by lack of concentrated recharge from overlying beds.
II.III. Coastal aquifer and density dependent flow model
When groundwater quality is contaminated, it is mainly due to the spread of a polluting solute through an aquifer. A known examples are seawater intrusion and waste disposal in deep salt formations. These two examples are well studied by (Buès and Oltean 2000; H. J. Diersch 1988; Frind 1982; G. Galeati et al. 1992; Herbert et al. 1988a; Ekkehard Holzbecher 1998; Huyakorn et al. 1987; Johns and Rivera 1996; Olaf Kolditz et al. 1998; L. F. Konikow et al. 1997; C. Langevin et al. 2005; Misut and Voss 2004; Oldenburg and Pruess 1995; Putti and Paniconi 1995; Segol and Pinder 1976; Segol et al. 1975a; Van Camp et al. 2014; Volker and Rushton 1982; C. I. Voss 1984a; C. I. Voss and Souza 1987; A. Younes et al. 1999). Freshwater stored in coastal aquifer is extensively extracted through pumping wells (Barlow and Reichard 2010; Post 2005) due to high water demand in coastal area (touristic, industrial and public use). In order to enhance freshwater security and to avoid contamination of our reserves (H. J. G. Diersch and Kolditz 2002; Post 2005), seawater intrusion becomes a topic of great interest for hydrogeologist. During the last few decades, hydrogeologist provided a deeper understanding of the prediction, processes, investigative tools and management of such systems (Thomas E. Reilly and Goodman 1985; Cheng and Ouazar 2003). During seawater intrusion seawater and freshwater mix more quickly, this phenomenon which called density effect, enhanced solute transport which can accelerate the hydrodynamic mixing of different fluids within shorter time scales and over longer distances than that caused by diffusion alone (H. J. G. Diersch and Kolditz 2002; C. Simmons 2005; C. T. Simmons et al. 2001). Density effect is the term that classifies the flow pattern influenced by density differences in the fluid system (Ekkehard Holzbecher 1998). Many works use the term density-driven flow since most of variable density flows in a porous medium in nature are generated by density differences. Variable density flow analysis is associated not only with the density itself but also with the parameters that affect or change the density of the fluid. These parameters can be temperature, concentration (salinity) or both simultaneously. Except density effect, other factors and processes associated with seawater intrusion need to be taken into account , including density variations induced by seawater penetrating into the freshwater (Abd-Elhamid and Javadi 2011; Thomas E. Reilly and Goodman 1985), porous media heterogeneity (H. J. G. Diersch and Kolditz 2002; Kerrou and Renard 2010; X. Li et al. 2009; Mulligan et al. 2007; Konz et al. 2009), the width of the mixing zone (Abarca et al. 2007a), seasonal fluctuations in influx of seawater and groundwater (Michael et al. 2005; Taniguchi et al. 2006), tidal effects (Ataie-Ashtiani et al. 2001; Boufadel et al. 2011; H. Li et al. 2008; Lu et al. 2009; Robinson et al. 2007), dense contaminants leaking to the Aquifer (C. Simmons 2005), geothermal convection (W. S. Moore 2010; Wilson 2005) and the chemistry of the fluid and aquifer (Boluda-Botella et al. 2008; Rezaei et al. 2005; Santoro 2010) .
II.III.1 The most known situation of seawater intrusion
Seawater Intrusion generally can be categorized into one or more of several types of intrusion: horizontal and upward movement of the interface, downward leakage of brackish or saltwater from surface water (such as in estuarine environments), or saltwater upconing beneath a well field.
II.III.1.1 Horizontal intrusion
As shown in Fig. II.6 (a) horizontal intrusion occurs as the saltwater from the coast slowly pushes the fresh inland groundwater landward and upward. This type of intrusion can be regional in scale, and results in the characteristic “wedge” of saltwater at the bottom of an aquifer. Its cause can be both natural (due to rising sea levels) and pumping of freshwater from coastal wells (Misut and Voss 2004; Van Camp et al. 2014). There is always an interface between the saltwater offshore and the freshwater onshore. This interface can sometimes be relatively sharp, with little or no transition or diffusion zone, this have been seen on Long Island (Misut and Voss 2004).
II.III.1.2 Downward leakage of saltwater from surface water
Pumping from coastal wells can also draw saltwater downward from surface sources such as tidal creeks, canals, and embayments. This type of intrusion is shown in Fig. II.6 (b), is usually more local in nature. It typically occurs within the zone of capture of pumping wells where significant drawdown of the water table causes induced surface infiltration. This type of intrusion has occurred in areas in southern Everglades of Florida (C. Langevin et al. 2005), where leakage between the surface and subsurface is locally important in the wetland.
Figure II.6 Seawater intrusion situations, downward leakage of saltwater from surface water (a), Induced downward movement of brackish surface water (b), saltwater upconing beneath a supply well (c) (Cheng and Ouazar 2003)
II.III.1.3 Upconing below pumping wells
Upconing (Fig. II.6 (c)) below pumping wells was studied numerically by (H. J. Diersch et al. 1984; Ekkehard Holzbecher 1998; T. E. Reilly and Goodman 1987; Johannsen et al. 2002). Saltwater upconing is described by (T. E. Reilly and Goodman 1987) as the movement of saltwater from a deeper saltwater zone into the fresh groundwater in response to pumping at a single well. This is generally a more local intrusion problem, experienced by individual wells or wellfields. It requires that the pumping well be screened in freshwater overlying saltwater. At a certain pumping rate, a stable cone in the interface can develop below the well screen, but will not rise to the well. At increased rates of pumping, the cone can become unstable and the interface will rise abruptly toward and into the well, causing the well discharge to become saline. It is a widespread problem reported in USA that is observed at Georgia, Michigan, Minnesota, Missouri, Nebraska, New Mexico, North Dakota, South Carolina, and Utah (Krause and Clarke 2004) as well as globally in other countries such as Germany (H. J. Diersch and Nillert. 1987; H. J. Diersch et al. 1984), the Netherlands (Huisman 1954), Israel (Dagan and Bear 1968; Schmork and Mercado 1969), Japan (Hosokawa et al. 1990), Sarir Tazerbo wellfield in Libya (Laabidi et al. 2010) .
II.III.2 Survey of Numerical model and benchmarking density dependent flow models:
Numerical simulations (Abd-Elhamid and Javadi 2011; Albets-Chico and Kassinos 2013; H. J. G. Diersch and Kolditz 2002; D. Henry et al. 2012; Thorne et al. 2006; Bouhlila 1996), field studies (Walther et al. 2012; Zimmermann et al. 2006) and laboratory studies (Jakovovic et al. 2011; Konz et al. 2009; Schincariol and Schwartz 1990) have shown the important effect of density and viscosity gradients caused by the variation in concentration on flow and solute transport in porous media. When there is a difference between the densities of two mixed waters, density-driven flow can manifest itself in the form of lobe-shaped instabilities or fingers; this is the main cause of the free convection (H. J. G. Diersch and Kolditz 2002; D. Henry et al. 2012; Olaf Kolditz et al. 1998; C. T. Simmons et al. 2001) .
The few analytical solutions for variable density flow problems create a significant dependence on numerical codes for simulating seawater intrusion. Numerical modeling is typically the tools to understand the influence of the processes listed above, and well predict the behavior of such systems. Such computer models have become essential tools to interpret real situation and laboratory studies of seawater intrusion (Adrian D. Werner et al. 2013). Most codes solve a coupled system of variable density flow and solute transport equations, which is numerically complicated and a computational effort is needed. The difficulty is mainly related to the solute transport part, which requires a fine discretization of the aquifer (C. I. Voss and Souza 1987; H. J. Diersch 1988; H. J. Diersch et al. 1984; H. J. G. Diersch and Kolditz 2002). Common techniques such as finite volume, finite element or finite difference methods perform well for the solution of the flow equations, but they result in numerical dispersion when applied to solve the solute transport equations (H. J. G. Diersch and Kolditz 2002). Among the large number of studies carried out during the last years, the following basic studies were highlighted:
(H. J. Diersch 1988; H. J. Diersch et al. 1984; H. J. G. Diersch and Kolditz 2002; Rifai et al. 1956) showed that the dispersion coefficient is proportional to the velocity and it is much greater in the direction of the velocity than in the lateral direction.
(H. R. Henry 1964a) developed the first analytical solution including the effect of dispersion in confined coastal aquifers under steady-state conditions. The author assumed a constant dispersion coefficient in the aquifer. Henry’s problem was studied by Lee and Cheng (1974) in terms of stream functions.
(Segol et al. 1975a) developed the first transient solution based on a velocity-dependant dispersion coefficient. They used Galerkine finite element method to solve the set of nonlinear partial differential equations describing the movement of a saltwater front in a confined coastal aquifer.
(Simpson and Clement 2003) Improving the worthiness of the Henry problem as a benchmark for density-dependent groundwater flow models. They present an alternative solution to the standard Henry saltwater intrusion problem where reduction of freshwater recharge is tested. The reduction of the freshwater recharge means that the density effects play a more significant role in determining the isochlor distribution (Kerrou and Renard 2010).
(Abarca et al. 2007a) developed a more generic description of the seawater intrusion problem exemplified by the Henry problem that includes both velocity dependent dispersion and anisotropy in hydraulic conductivity.
(Pool and Carrera 2010) showed the significant contribution of transverse dispersion in creating mixing in the interface.
We can now find several models dealing with the problems listed above with broadly similar methodologies. The best known and the most widely used by the hydrogeologists community are as follow:
SUTRA is a well-documented 2D finite element code (C. I. Voss 1984a) for simulating variable-density groundwater flow model throughout the world by (Souza and Voss 1987; C. I. Voss and Souza 1987). It can simulate density dependent groundwater flow with energy transport (thermal) or chemically reactive (single specie) solute transport.
HST3D is a 3D finite difference code that can simulate heat and solute transport (Kipp 1986). Modeling a large-scale coastal hydrogeologic system with a non-uniform density distribution with HST3D appears to be rather complicated when small longitudinal dispersivities should be applied. For heat transport applications, this code seems to be much more appropriate; for instance, it often used for modeling geothermal systems in the Netherlands.
SWICHA is a 3D finite element code (Huyakorn et al. 1987), it can simulate variable-density fluid flow and solute transport processes in saturated porous media. The applications range from simple 1D to complex 3D, coupled flow and solute transport. The solute transport equation may not converge to a solution if the Peclet number is exceeded in an element. To solve this problem, SWICHA offers a trick at the user's option: artificial dispersion is added to the solute transport equation matrix. Spatial oscillations are then suppressed and the critical Peclet number is no longer exceeded in that element.
METROPOL (MEthod for the TRansport Of POLlutants) is a 3D finite element code for simulating groundwater flow with varying density and simultaneous transport of contaminants (Sauter and Leijnse 1993). METROPOL is developed by the Dutch National Institute of Public Health and Environmental Protection RIVM. It has been applied to simulate safety assessments of the geological disposal of radio nuclear waste in salt formations.
FAST-C(2D/3D) (E. Holzbecher 1998) can be applied to model density driven flow in porous media using the stream function formulation. Density variations are caused by temperature or salinity gradients. Steady state as well as transient modelling of confined aquifers are possible.
FEFLOW is a 3D finite element computer code, which employs the finite element method (H.J.G. Diersch 1996). The governing partial differential equations describe groundwater flow, where differences in density affect the fluid flow is implemented. Fluid density effects are caused by contaminant mass as well as temperature differences simultaneously, inducing thermohaline flow.
D3F (Distributed Density Driven Flow) is one of the sophisticated, computer codes to simulate density-driven groundwater flow (Fein 1998). It is based on the volume element method. Permeabilities, porosities, values for boundary and initial conditions can be defined as constants as well as simple spatial user-defined functions. In addition, a stochastic permeability distribution technique is available to consider small-scale heterogeneities. Fluid density and viscosity are functions of salt concentration and temperature. The code can run on serial as well as massively parallel computers.
MOCDENS3D (G. H. P. Oude Essink 2001) is MOC3D (L.F. Konikow and Goode 1996), but adapted for density differences. Buoyancy terms are introduced in the vertical effective velocities to take into account density effects. The adapted MODFLOW module solves the equation of density dependent groundwater flow. Solute transport is modeled by the method of characteristics (MOC) module in MOCDENS3D. The particle tracking method is used to solve advective solute transport and the finite difference method is used for hydrodynamic dispersive transport.
SEAWAT (Guo and Langevin 2002a) is a combination of MODFLOW and MT3DMS (Zheng and Wang 1999). It is designed to simulate 3D variable-density ground water flow and solute transport. The program was developed by modifying MODFLOW subroutines to solve a variable-density form of the ground water flow equation and by combining MODFLOW and MT3DMS into a single program.
Besides the programs listed above, as presented in Table II.2, we can find in the work presented by (Adrian D. Werner et al. 2013) a summary of the popular seawater intrusion codes, their basic features, and model testing.
Table II.2 Popular density dependent flow codes (Adrian D. Werner et al. 2013)
FD – finite difference; FE – finite elements; IF – interface flow; S – saturated flow only; SU – saturated–unsaturated flow; GUI – dedicated graphical user interface available. E – Elder problem; H – Henry problem; C – HYDROCOIN salt-dome problem; R – rotating fluids problem.
Density-dependent flow models must be verified in a particular series of benchmarks. Standard tests are Henry, Elder, salt dome, and the convection cell problems. (H. J. G. Diersch and Kolditz 2002), (C. Voss et al. 2010) provide a discussion on the applicability of different tests to variable density flow codes. The first benchmark used was presented by Henry (1964) who solved the governing equations of the density-dependent flow in a non-dimensional form and presented a semi-analytical solution for the steady state distribution of salt concentration and stream functions. Although it does not exhibit significant density gradients, the Henry problem has been widely used as a benchmark of density-dependent groundwater flow (Abarca et al. 2007a; H. J. G. Diersch and Kolditz 2002; Simpson and Clement 2003; C. I. Voss 1984a). Other popular benchmark problems include the HYDROCOIN salt dome problem (L. F. Konikow et al. 1997), salt lake problem (C. T. Simmons et al. 1999), salt pool problem (Johannsen et al. 2002; Oswald and Kinzelbach 2004), rotating fluids problem (Bakker et al. 2004), and the Henry and Hilleke problem (C. D. Langevin et al. 2010).
II.IV Geochemical Settings
Carbonates are minerals characterized by the ion. The most abundant ones are calcite (CaCO3), aragonite (CaCO3), and dolomite (Ca,Mg)(CO3)2. Dolomite is a frequent mineral in various sedimentary rocks such as dolomitic limestones, evaporates, etc. as in the metamorphic rocks, which are derived from them (Morel, 1994). The carbonated rocks are sedimentary ones formed at least by 50% of carbonate (calcite, dolomite, and aragonite).
Dissolution or precipitation of carbonates plays an important role in nature, such as in the chemistry of seawater, in the geochemical evolution of freshwater and the evolution of karst environment. Many investigations have been carried out in last few decades in order to provide a deeper understand of these process (reference). Dissolution and precipitation reactions can be described by the mass-action law as reversible, heterogeneous reactions. In general, the solubility of a mineral is defined as the mass of mineral, which can be dissolved or precipitated within a standard volume of the solvent. As for any mineral, chemical equilibrium reaction is described by the solubility product constant K and the ion activity product IAP.
Dissolution and precipitation of calcium carbonate are determined by the chemistry of the system H2O-CO2-CaCO3. The chemical reaction is usually given by the following equations:
K1=10-3.5
K2 =10-10.25
KH =10-2.8
Kw =10-14
Kc =10-8.4
Equation (R1) and (R2) represent the ionizations of carbonic acid, Equation (R3) represents the statement of Henry’s law, equation (R4) describes the dissociation of water, equation (R5) describes the solubility product for calcium carbonate and finally equation (R6) represents the electroneutrality condition.
Whether a water is undersaturated, saturated, or supersaturated with calcite can be determined by comparing the equilibrium constant K of the mineral of interest to the ion activity product (IAP) of that mineral in the water. The IAP of calcite is computed as follows: and represents the saturation state which is the ratio of IAP to the solubility product K and is defined as the saturation index SI of calcite which represents the equilibrium state of a solubility reaction.
If SI<0, then the solution is undersaturated with respect to calcite and dissolution occurs.
If SI=0, then the solution is in equilibrium
If SI>0, then the solution is supersaturated with respect to calcite and precipitation occurs.
II.IV.1 Kinetic reactions in sedimentary basin
Several semi-empirical solutions are available in literature describing kinetics of calcite or dolomite dissolution and precipitation. But in most cases, they are based on kinetic constants determined in laboratory in controlled conditions for a simple system (C.A.J. Appelo et al. 1984; Arvidson et al. 2003; Buhmann and Dreybrodt 1985; Busenberg and Plummer 1982; Dreybrodt 1980; Herman and White 1985; Morse and Arvidson 2002; L. N. Plummer et al. 1979)
Sedimentary basins can contain close to 20% by volume of pore fluids commonly classified as brines. The water in a basin can become undersaturated with respect to calcite as a result of convection, dispersivity, or injection of CO2 in geological reservoir.
The chemical reactions between a solid and a solution are quantified by the mass change of the solid and the species concentrations in the solution. Many authors have been working toward an understanding the kinetic of calcite dissolution in low temperature aqueous systems. The most commonly used equation describing the kinetic of carbonate mineral dissolution and precipitation was presented by (Morse and Berner 1972) as follows:
, (II.1)
where is the rate in (µmol m-2 h-1 ), m is moles of calcite, t is the time, A is the total surface area of the solid, V is the volume of the solution, is the rate constant and n is a positive constant describing the order of the reaction.
In (Morse and Arvidson 2002) there is a detailed presentation of a modified general form of the equation providing a bridge between thermodynamics and Kinetic as follows:
, (II.2)
where A is the reactive surface area, Ea is the apparent activation energy, T is the absolute temperature, R is the ideal gas constant, and ai is the activity of the ions in the solution.
Many investigations have been presented in the literature to analyze the effect of the thermodynamic and chemical parameters in the rate. We can find in (Gledhill and Morse 2006) a complete analysis of the multiple regression (MR) model of the form as follows:
, (II.3)
Where T is the temperature in °C, pCO2 is given in bars, I is the ionic strength and ai are activities of calcium and magnesium. According to (L. N. Plummer et al. 1979) and in a simple chemical system, the dissolution or the precipitation rate of calcite is a function of the activities of the species in the solution, temperature, and pCO2 and is described by the following equation:
. (II.4)
One of the important factor is the pH of the solution, many investigations were carried out during the last few years to understand and evaluate the rate of the dissolution and precipitation of a solid during the chemical reaction, we summarize in Fig. II.7 some of the previous works that focus on the evaluation of the pH effect on the calcite dissolution rate.
Figure II.7 Effect of pH on the calcite dissolution rates from previous studies
(Morse and Arvidson 2002).
(Morse and Arvidson 2002) present in their work a series of chemical and physical processes of the kinetic chemical reaction that can be written in seven steps:
Step1: diffusion of the reactants to the solid;
Step2: adsorption of the reactants on the solid surface;
Step3: migration of the reactants to a reactive area in the solid surface;
Step4: chemical reaction between the reactants and the solid surface;
Step5: migration of the products (results of the chemical reaction) away from the reaction area;
Step6: desorption of the products to the solution;
Step7: diffusion of the products.
II.IV.2 Mixing of carbonate waters
The theory of mineral saturation by mixing has been extensively used to interpret dissolution or precipitation process (Carrera et al. 2004; Chu et al. 2005; Cirpka et al. 1999; Cirpka and Valocchi 2007; Corbella et al. 2003; De Simoni et al. 2007; Emmanuel and Berkowitz 2005; Hanor 1978; Knutson et al. 2005; Rezaei et al. 2005; Romanov and Dreybrodt 2006; Rueedi et al. 2005; Sanford and Konikow 1989a; Stoessell et al. 1989). Alterations in mineral composition or precipitation or dissolution can occur during mixing of different waters. One common example is the mixing of saltwater and freshwater of coastal carbonate aquifer. (L. N. Plummer 1975a) and (Wigley and Plummer 1976a) have provided a clear analysis of geochemical process during mixing of freshwater and seawater in carbonate coastal aquifer. The most important result that they reached is that mixing of freshwater saturated with respect to calcite with seawater that is supersaturated with respect to calcite can lead to a solution undersaturated with respect to calcite. The resultant mixture is strongly dependent on the algebraic effect, on the temperature, on the nonlinear dependence of activity coefficients on ionic strength, on the nonlinear variation of equilibrium constant with temperature, on the ionic strength effect, and most importantly on the partial pressure of CO2 (pCO2). Similarly, Bogli (1964) showed that mixing of two carbonate groundwaters which are both saturated with respect to calcite can lead to undersaturated conditions if the original groundwaters have different pCO2.
The algebraic effect
Mixing of two fluids saturated with calcite but with different Ca concentrations under closed system conditions typically produces a fluid that is supersaturated with calcite. This occurs because fluids equilibrated with calcite that have different Ca concentrations also have different carbonate activities (Corbella et al. 2003). In general if ion activity coefficients are assumed constant, the ion activity product IAP of any mixture is always higher than the equilibrium constant. This causes supersaturation of mixtures, even if the end members are in equilibrium (Fig. II.8 (A)).
Salinity and ionic strength effects
As shown in Fig. II.8 (B) ionic strength effects occurs when two solutions with different salinity are mixed together. Mixing of two fluids saturated with calcite and chemically similar except for their salinities results in a solution that is undersaturated in calcite and may lead to dissolution of calcite (Corbella et al. 2003). This occurs because calcium activity is lower in solutions of intermediate salinity than in solutions with a very high or low salinity (Wigley and Plummer 1976a).
Figure II.8 Schemes illustrating three of the most relevant factors affecting calcite saturation of a mixture of two calcite-saturated solutions. A) the algebraic effect. B) the salinity effect (the labels % indicate the proportion of of a 1m NaCl-solution in the mixture). C) pH and pCO2 effect.
Effect of pH and pCO2 (Fig. II.8 (C))
The saturation of a mixture also depends on the pH and pCO2 because of the strong dependence of the carbonate species on pH (Fig. 1C). Thus, the CO3 concentration will be higher under slightly basic conditions than under acidic ones. Therefore, the mixture of slightly acidic fluids will be more easily subsaturated than the mixture of more basic ones.
Effect of temperature
Mixing of fluids saturated with calcite but with different temperatures may result in a solution that is undersaturated or supersaturated with respect to calcite depending on the chemistry of the end members (Corbella et al. 2003; Wigley and Plummer 1976a). The temperature effect is relatively minor except where there is a significant temperature difference between the fluids such as occurs when a hot geothermal fluid mixes with shallow groundwater or seawater, or a deeper basinal fluid is able to migrate relatively quickly to a shallower unit.
Many models have been published in the last few decades that describe the geochemistry of natural water: TEQUIL (Greenberg and Møller 1989), GIMRT (C.I. Steefel and Yabusaki 1996), PHREEQC (Parkhurst and Apello 1999a), EQL/EVP (Risacher and Clement 2001), EQ3/6 (Wolery 1992), CHESS (J. Van Der Lee 1998), Geochemist's Workbench (Bethke 1994), and GEMS3K (Kulik et al. 2013). In fact, the total number of programs that are available both commercially and in public domain of geochemical processes in subsurface systems is likely to be significantly larger than the number of programs listed above.
II.V Reactive transport modeling in variable density flow
An important aspect of reactive transport modeling is the appropriate evaluation of activity coefficients both under dilute and highly saline conditions. Most reactive transport codes are based on the Davis equation or some form of the extended Debye-Huckel equation to determine the activity of dissolved ions as a function of ionic strength. However, these activity relationships are only applicable in relatively dilute solutions up to an ionic strength little to that of seawater. Highly saline solutions and brines found in sedimentary basins require an alternative approach for calculating the activities of dissolved ions. One of the most common approaches is the Pitzer ion interaction model (Pitzer 1973; Pitzer et al. 1984). The Pitzer ion interaction model is based on a virial expansion, which reduces to a modified Debye-Huckel formulation at low concentrations (Peter C. Lichtner and Felmy 2003). In addition to the effects of the formation of complexes and ion pairing, activity coefficients calculated using the Pitzer formalism account for electrostatic interactions and ion hydration effects (Peter C. Lichtner and Felmy 2003). However, considering the significant salinity gradients that may be present, a coupled treatment of geochemical reactions and density dependent flow and transport is needed and fewer codes meet these requirements (V. Freedman and Ibaraki 2002; Henderson et al. 2009; Kim et al. 2004; Peter C. Lichtner et al. 2004; Mao et al. 2006b). Although these codes are theoretically suitable, they have not been used to date to evaluate the long-term geochemical evolution in sedimentary basins. The effect of dissolution-precipitation reactions on porosity and permeability has also been implemented in reactive transport codes (Henderson et al. 2009; Peter C. Lichtner et al. 2004; C. I. Steefel and Lasaga 1994). In addition, the capability of these models to deal with strong salinity contrasts, locally diffusion controlled transport, and relatively complex geometries needs to be evaluated.
There is a large number of studies where the chemistry is highly complex, with a large number of mobile species involved in the reactions (Bolton et al. 1996; V. Freedman and Ibaraki 2002; V. L. Freedman and Ibaraki 2003; Mayer et al. 2002; Saaltink et al. 2001; Sanford and Konikow 1989a; D. Schäfer et al. 1998; W. Schäfer and Therrien 1995; Spycher et al. 2003; C. I. Steefel and Lasaga 1994; C. I. Steefel and Lichtner 1998; Walter et al. 1994; Yeh and Tripathi 1989; Zysset et al. 1994). In recent decades, reactive transport models become important tools for enhancing the understanding of the hydrogeochemical process in groundwater systems. There are several reactive transport models available in the literature that couple flow and transport code with geochemical one. One of the biggest issues in numerical implementation of such problem is the scheme to couple the transport process and the geochemical reactions. The two commonly methods used for coupling transport and reaction process in the literature are the one step approach (also called the global implicit approach or the direct substitution approach:DSA), and the two-step approach (also can be called sequential approach or operator-splitting approach).
Global implicit method, One-step or Direct Substitution Approach (DSA): in this method the governing equations including both transport and reaction terms are solved simultaneously which leads to a fully coupled system of equations. The chemical equations are substituted into the transport equations and it is solved by applying the Newton-Raphson method. This method has been used by (V. Freedman and Ibaraki 2002; Glassley et al. 2003; Grindrod and Takase 1996; Mayer et al. 2002; Miller and Benson 1983; M. W. Saaltink et al. 2004; C. I. Steefel and Lasaga 1994; Tianfu Xu et al. 2006b).
Operator splitting technique, where the calculation for transport and reaction equations is decoupled and the partial differential equation for chemical transport and chemical reactions are solved sequentially. This method has been used by (C. A. J. Appelo and Willemsen 1987; Centler et al. 2010; Engesgaard and Kipp 1992; Engesgaard and Traberg 1996; Kosakowski and Watanabe 2014a; Le Gallo et al. 1998; Liu and Narasimhan 1989; Mao et al. 2006b; Parkhurst and Apello 1999a; Parkhurst et al. 2004; Poulet et al. 2012; Prommer 2002; J. Samper et al. 2000; Javier Samper et al. 2009; D. Schäfer et al. 1998; Šimůnek and Suarez 1994; C.I. Steefel and Yabusaki 1996; Tebes-Stevens et al. 1998; Jan van der Lee et al. 2003; Walter et al. 1994; Wernberg 1998; Tianfu Xu et al. 1999; Yeh and Tripathi 1990). The sequential approaches including SIA and SNIA are among the proposed approaches in this category. In the case of sequential iterative approach (SIA), the transport equations are solved first, then using the results obtained from the transport step, values of chemical concentrations are corrected (MacQuarrie and Mayer 2005). This procedure is repeated for every iteration within a time step. Sequential non-iterative approach (SNIA) is a special case of SIA where the chemical equations are solved only once per time-step, i.e., after the transport step has achieved converged results (Walter et al. 1994).
Several authors have discussed the advantages and disadvantages of the various solution approaches used in reactive transport modeling (Mangold and Tsang 1991; Saaltink et al. 2001; C. I. Steefel and Lasaga 1994; Yeh and Tripathi 1989). An advantage of the two-step method is a higher flexibility when dealing with a complex chemical system, which makes it easy to add or delete a process.
(Carl I. Steefel and MacQuarrie 1996) compared the results of simulations in simple cases using different sequential coupling methods. They reported that SIA sometimes gives the smallest error to CPU time ratio, although in other cases, SNIA is more efficient. However, the accuracy of SNIA is highly dependent on space and time discretisation and the type of chemical reactions. (Tianfu Xu et al. 1999) also employed SNIA in addition to SIA method and compared the accuracy and numerical performance of using several test cases. They demonstrated for a specific case that SNIA is more efficient and requires less CPU time. However, numerical dispersion was noticed and the accuracy also reduced. Therefore they concluded that there is no clear proof that SNIA is better than SIA, and it mainly depends on the choice of the chemical reactions and also on the time and space discretisation. (Javier Samper et al. 2009) showed that numerically SNIA is two to three times more efficient than SIA. Generally, the sequential approach is easier to program and the memory requirements are less than those for a fully coupled method. Advantages of the one-step method are the simultaneous treatment of all processes in time, physical space, and chemical reaction space. As well, the global convergence properties of the one-step method may be better than for the two-step method, and it is sometimes possible to take larger time steps compared with those in the two-step approach.
A number of models and codes have been developed during the last few decades (table II.3). Many difference between these models can be identified such as the coupling method, numerical method (finite element, finite difference or finite volume method), dimension of the model (1-D, 2-D, or 3-D), the capability to simulate density driven, protocol and model used to simulate permeability-porosity feedback (permeability-porosity relationship into reactive model), capability to simulate saturated/unsaturated transport, and method to calculate ion activity and ion speciation (Davies model, Debye-HucKel model, and Pitzer activity model). In this regard, table II.3 presents a summary of the coupling schemes implemented for solving the coupled transport and geochemical reactions in popular reactive transport codes, and in the following we present a brief description of their capabilities and features:
HYDROGEOCHEM (Yeh and Tripathi 1990) is one of the first comprehensive simulators of subsurface reactive transport. HYDROGEOCHEM is a 2-D numerical model of coupled flow, thermal, transport and geochemical reaction in variable saturated porous media. The mathematical model was solved using the finite element method. The coupling method between the transport and the geochemical reaction was maintained using the sequential iterative/non-iterative approach. It covers all geochemical reaction such as aqueous complexation, adsorption, precipitation-dissolution and ion-exchange.
The ARASE model (Grindrod and Takase 1996), a 2-D numerical model for coupled flow, thermal, transport and geochemical reaction in variable saturated porous media. It was one of the first models to fully couple flow, transport and geochemical reactions, and the porosity-permeability relationship (permeability feedback created by volumetric minerals change during dissolution-precipitation reactions). The mathematical model was solved using the global implicit approach.
The RETRASO code (Reactive T RAnsport of SOlutes) was developed by the Technical University of Catalania (M. W. Saaltink et al. 2004; Saaltink et al. 1998) in order to model equilibrium and kinetic reactions for non-conservative dissolved species. This was originally a 2-D finite element code for saturated porous media. RETRASO is based on the global implicit approach to solve the coupled transport process and geochemical reactions. The result is carried out using a Newton-Raphson method. The new version of this code was presented by (M. W. Saaltink et al. 2004). The model is the coupling of the old code with the multiphase flow and heat code called CODE-BRIGHT model (Olivella et al. 1996). The new code can handled both saturated and unsaturated flow, heat transport and reactive transport in both liquid and gas phases.
(Le Gallo et al. 1998) developed DIAPHORE model, a 3-D model to simulate water-rock interaction and porosity-permeability evolution during geochemical reactions such as precipitation-dissolution reactions. The model was used to simulate porosity-permeability changes in an oilfield. The authors were able to fit the observed values for concentrations and permeability in the modeled oilfield. The mathematical problem was solved using a finite volume method and the coupling scheme between the transport process and the geochemical reactions was implemented using a sequential iterative approach.
(Prommer 2002) presented a computer code PHT3D for reactive transport in porous media, coupling MODFLOW/MT3DMS for transport and PHREEQC for chemical reactions. The model uses a sequential non-iterative approach to couple the transport and the chemical reactions. The concentration distribution results of the transport problem of the first time step were used as an input file for PHREEQC for calculating chemical reactions such as ion-exchange, dissolution-precipitation… This protocol is used for the next time steps until reaching the simulation time.
(V. Freedman and Ibaraki 2002) used the global implicit approach to develop a reactive transport model with variable density flow called DART. The model was developed by combining the 3D density dependent flow and transport code MITSU3D (Ibaraki 1998), with the geochemical code GIMRT (C.I. Steefel and Yabusaki 1996). The model was used to calculate the effect of the density dependent flow and geochemical reactions in porous media.
COMPASS model was presented by (Seetharam 2003), the geochemical reactions were calculated using MINTEQA2 model (Allison et al. 1991). Both sequential iterative and non-iterative approach have been used in this model.
(Mao et al. 2006b) presented a 2-D finite difference model called PHWAT. PHWAT is a combination between the density dependent groundwater flow and solute transport model, SEAWAT (Guo and Langevin 2002a), and the geochemical reactions model, PHREEQC-2. The coupling between the transport process and the geochemical reactions was implemented using a sequential iterative/non-iterative approach. The model was evaluated using experimental data of seawater intrusion in a column filled with sandy sediment. The model was shown a capability to model physically unstable flow and that numerical instabilities were suppressed (Mao et al. 2006b).
(Centler et al. 2010) present GEOSYSBRNS, a reactive transport code by combining the BRNS geochemical code (Aguilera et al. 2005) to the flow and transport simulator GEOSYS (Olaf Kolditz and Bauer 2004). BRNS can handle kinetically the biogeochemical reactions, and GEOSYS has the capability to simulate multidimensional finite element subsurface flow and transport. The coupling scheme is based on a sequential non-iterative approach.
ESCRIPT-RT a reactive transport code for fully saturated porous media developed by (Poulet et al. 2012) which is based on a finite element method coupled to a Gibbs Energy Minimization solver, to an equation of state to calculate fluid proprieties, and to thermodynamic database to calculate rocks proprieties. This code uses a sequential scheme to solve temperature, pressure, mass transport and chemical equilibrium.
(Kosakowski and Watanabe 2014a) developed OpenGeoSys-Gem, a reactive transport code used to solve chemical equilibria. The model is developed by combining the flow and transport model OpenGeoSys (O. Kolditz et al. 2012) and the GEM based chemical solver GEMS3K (Kulik et al. 2013). The coupling scheme is based on a sequential non-iterative approach.
Table II.3 Popular reactive transport codes
-SIA : Sequential iterative approach -SNIA: Sequential non-iterative approach -DSA: Direct substitution approach -FE: Finite element -FD: Finite difference -FV: Finite volume
II.VI Reactive transport modeling in carbonate coastal aquifer
This section aims to give a short overview of related studies dealing with reactive transport in carbonate coastal aquifers and CO2 storage associated with geochemical reactions, the most known works are:
(Sanford and Konikow 1989a) provided a quantitative estimate of the amount of calcite dissolution expected in mixing zone under typical hydrodynamic and geochemical conditions. To assess the effects of calcite dissolution on the rock matrix, the authors proposed a procedure to quantify the effect on the rock proprieties (porosity and permeability). They used a two-step method consisting basically of (1) transport of TDS and dissolved calcite, (2) evaluating at each node the dissolution potential (derived from the mixin ratio), and (3) dissolving the dificit. They used the geochemical reaction model Phreeqc (Parkhurst and Apello 1999a) to calculate the amount of dissolved calcite in the freshwater seawtaer mixing zone. After calculating the porosity changes they show that the increases of porosity due to the quantity of dissolved calcite released from the pore volume could lead an increase of permeability which affects the flow field and solute transport.
(Rezaei et al. 2005) extended the study of (Sanford and Konikow 1989a) simulated the porosity development in coastal aquifer due to calcite dissolution. They found that maximum dissolution occurs neither where transport enhance mixing nor where the maximum undersaturation. The authors used RETRASO code (M. W. Saaltink et al. 2004) which solve simultanuosly the aqueous complexation reaction, the dissolution and precipitation reactions and the flow and transport problem using a global implicit method. They found through the simulations two locations of maximum dissolution as obtained by (Sanford and Konikow 1989a), one near the aquifer bottom and one near the discharge boundary. The authors did not take into account in their simulations the feedback effect of porosity and permeability which would enhance further mixing.
(Romanov and Dreybrodt 2006) used the analytical solution of (Phillips 1991)to set up an alternative modeling approach to determinate the porosity change in the saltwater-freshwater mixing zone of coastal carbonate aquifers. As a first step, they solve the advection transport equation for salinity s using the SEAWAT code (Guo and Langevin 2002a) and the second step consists in calculating the dissolution capacity in the mixing zone by estimating the calcium equilibrium concentration that is obtained as a function of salinity using PHREEQC (Parkhurst and Apello 1999a). They applied this approach for a stationary salinity distribution. Nevertheless porosity change caused a change of permeability which causes variations of salinity distribution. In this regard such processes should be treated step by step for a non-stationary distribution of salinity.
(De Simoni et al. 2005a) presented a new procedure (general form of Phillips’s concept) to solve multicomponent reactive transport problems in porous media where homogenous or classical heterogeneous reactions occur (Rubin 1983a). The methodology consists of a protocol composed of four steps: (1) definition of mobile conservative components, (2) solving the transport equation for conservative species, (3) evaluation of the concentration of aqueous species and (4) substituting the latter into the transport equation to evaluate the reaction rate. The reaction rate is proportional to mixing rate where u the vector of conservative components concentration and D is the dispersion tensor. They used this approach to solve the pseudo-stationary mixing zone presented by (Sanford and Konikow 1989a). Their results is in good agreement with the result produced by the numerical code RETRASO (M. W. Saaltink et al. 2004).
(Singurindy and Berkowitz 2003a) performed a laboratory experiment to understand the evolution of hydraulic conductivity and flow patterns caused by dissolution of dolomite and precipitation of calcium carbonate in porous dolomite rock by using a column flow experiments. They conclude that long term experiments show that the evolution of hydraulic conductivity in carbonate rock, due to competition among flow, precipitation and dissolution processes, can display a wide variety of possible behaviors. In particular, a range of injected H/SO4 ratios and flow rates exists in which the effective hydraulic conductivity of the evolving carbonate rock samples displays oscillatory patterns. They conclude that higher flow rates cause a more rapid dissolution of the porous medium; in such cases, with dissolution dominating, highly conductive flow wormholes will develop. At slower flow rates, no wormhole formation will occur, but the porosity will vary in different parts of the domain. The authors propose the conclusion of their experiments as tools to provide information on the flow history of such formation.
(Singurindy et al. 2004a) performed a laboratory experiments to measure the total change of porosity in a 2-D flow cell filled with calcium carbonate particles (porosity of 33%). The dimension of flow cell is 0.16×0.16×0.01. The flow cell is connected to two inlet reservoirs. At the upside boundary, water in equilibrium with calcite at a carbon dioxide pressure of 1atm is injected with a total flux of 0.7×10-3 m3/day. At the left hand side and under the same flow boundary conditions, a solution of water containing NaCl is injected with concentration of 40g/l and of 80 g/l equilibrated with calcite and at pCO2 of 1 atm and 10-3.5 atm for the dissolution and the precipitation experiments respectively. They used a simple approach to verify the measured amount of dissolved or precipitated calcite based on Phillip’s approach (Phillips 1991) and they found a good agreement.
(V. Freedman and Ibaraki 2002) used a fully coupled density dependent flow and reactive mass transport model entitled DART (Ibaraki 1998). The model includes equilibrium reactions of aqueous species, kinetic reactions between the solid matrix and aqueous phases. A full coupling procedure of porosity and permeability change was implemented by the authors to quantify the impact of dissolution or precipitated mineral on the pore geometry. They study the effects of variable density fluid flow and reactive mass transport and their influence on generating instabilities in saturated porous media. They found that when permeability changes during dissolution reactions, the feedback loop does not have a significant impact on the plume shape.
(Katz et al. 2011) conducted a laboratory experiment in order to provide a quantitative analysis of calcite evolution patterns during mixing of two solutions of Na2CO3 and CaCl2 equilibrated with calcite. A 2D flow cell with dimension of 25cmx10cmx0.8cm in the x-, y- and z-direction respectively is used. The authors found that calcium carbonate precipitated along the mixing zone between the two solutions which flowed in parallel. The precipitation of calcium carbonate causes a decrease of porosity as pores in the mixing zone are filled by the precipitated mineral. In order to evaluate the effect of porosity change on the flow and transport problem, they measured the Ca2+ concentration in some selected sampling ports and used the RETRASO code (M. W. Saaltink et al. 2004) to numerically reproduce the results found in laboratory experiment. They found that the model appear to capture the main processes of the evolution of calcium concentration in flow domain.
(Audigane et al. 2007) modeled CO2 injection into the Utsira-Formation in Sleipner (North Sea, Norway) with geochemical reactions considering all important mineral phases. A 2D model was set up for an injection scenario of 25 years followed by a 10,000 years storage period. In this model the liquid concentration of Ca2+ and are assumed to be at equilibrium with a calcite mineral phase. The aim of this study was to quantify the mineralization process and therefore the overall fate of the injected CO2. The numerical simulator used is TOUGHREACT (Tianfu Xu et al. 2006b) which is a sequential iterative two phase flow, multi-component reactive transport code.
(T. Xu and Li 2013) created a realistic 1D geochemical flow model for the Frio-I Brine Pilot in Texas. Both, short and long-term effects of geochemistry were evaluated with TOUGHREACT (Tianfu Xu et al. 2006b) in this study, along with a model for the iron release due to the pH drop. The simulated results were compared to measurements taken at the site.
(Flukiger and Bernard 2009) set up a 3D fully-coupled reactive transport model with Stokes equations for flow description. The aim of the model is to interpret reactive percolation tests of carbonate samples on a laboratory scale.
III.I. Density-dependent flow and transport model
The equations governing the movement of a fluid through saturated porous medium subject to variable density conditions can be obtained from the mass and the momentum conservation principles. The following section presents briefly these governing equations (Bear 1988; H. J. G. Diersch and Kolditz 2002). For this study only two-dimensional flow is considered.
The first is the fluid mass conservation equation:
(III.1)
where kij[L2] is the intrinsic permeability, [L3L-3] is the porosity. [ML-3] and [ML-1T-1] are the fluid density and viscosity, respectively. These last two are dependent on the solute concentrations. p[ML-1T-2] is the fluid pressure, [ML-3T-1] is the source/sink term, mass flux per unit of volume, gj[LT-2] is the gravitational acceleration, and t[T] is time.
The GEODENS code describes transient convective and diffusive-dispersive transport of specie I, whose concentration is CI [MM-1], by the following equation:
(III.2)
where ui[LT-1] is the effective velocity in the porous media and is related to the Darcy velocity qi[LT-1] by the following expression:
(III.3)
and Dij[L2T-1] is the diffusion-dispersion term and is expressed as follows:
(III.4)
where is the Kronecker symbol, and are the longitudinal and the transversal dispersivity respectively, d0[L2] is the chemical diffusion, ui[LT-1] is the velocity vector component in the i direction, uj[LT-1] is the velocity vector component in the j direction, and u[LT-1] is the norm of the velocity vector.
The right-hand term of Eq.(III.1) can be formulated using the storage coefficient Ss[LT-2M-1] (C. I. Voss 1984a):
(III.5)
where β[LT2M-1] is the fluid compressibility and α[LT2M-1] is the solid matrix compressibility.
(III.6)
this can be formulated as follows:
(III.7)
by neglecting changes of as a function of C the following is obtained:
(III.8)
If Eq.(III.1) is combined with Eq.(III.2) by taking into account Eq.(III.3), this gives the non-conservative form of the transport equation:
(III.9)
Equations from (1) to (4) are completed by the state equations which describe the fluid properties:
(III.10)
The boundary conditions can be expressed as follows:
For the flow: In Г1 the pressure (p=p1) is imposed and in Г2 a flux (q=q1) is imposed.
For the transport: In Г3 the concentration (C=C1) is imposed and in Г4 a concentration total, convective or dispersive-convective flux (qc=qc1) is imposed.
And the following relationship is verified:
(III.11)
III.II. Geochemical Model
The main purpose of geochemistry model is to predict the minerals solubility, to quantify mineral phase variations in contact with an aqueous solution and to specify solute concentrations in this solution according to the thermodynamic conditions of the system (W. Stumm and Morgan 1970). Many models have been published in the last few decades that describe the geochemistry of natural water. The most common one is Phreeqc program (Parkhurst and Apello 1999a). Phreeqc is based on an ion-association aqueous model and has capabilities for speciation, saturation index calculations, batch-reaction and one-dimensional (1D) transport calculations involving reversible reactions. MINTEQ (Felmy et al. 1984) is a geochemical program to modelling aqueous solutions and the interactions of aqueous solutions with hypothesized assemblages of solid phases. MINTEQ applies the fundamental principles of thermodynamics to solve geochemical equilibria from a set of mass-balance equations. (Risacher and Clement 2001) published a FORTRAN program called EQL/EVP that simulates the evaporation of dilute waters as well as concentrated brines.
The prediction of mineral solubility is a complex problem in the geochemical modeling where two different methods can be applied. The first one is presented by (Harvie and Weare 1980) in which the mineral solubility in the brine is calculated by the minimization of Gibbs’ free energy under mass conservation and electrical neutrality conditions of the solution. The second one is based on the kinetic approach where the chemical reactions between a solid and a solution are quantified by the mass change of the solid and the species concentrations in the solution. Many authors have been working toward an understanding of the kinetic of mineral dissolution-precipitation in low-temperature of an aqueous systems (Gledhill and Morse 2006; Morse and Arvidson 2002; Morse and Berner 1972). The equation describing the rate is a function of the saturation index (Defined as the ratio of the ionic activity product IAP by the thermodynamic solubility product K) and is ranging from 0 to 1 for an unsaturated solution and takes a value greater than 1 for a supersaturated solution.
Highly saline solutions and brines found in sedimentary basins require a different approach for calculating the activities of dissolved ions. One of the most common approaches is the Pitzer ion interaction model (Pitzer 1973; Pitzer et al. 1984). The Pitzer ion interaction model is based on a virial expansion that takes into accounts various interactions in a solution such as cation-anion, cation-cation, anion-anion, cation-cation-anion, anion-anion-cation…
(Pitzer 1991) proposed an expression based on differentiating the excess Gibbs-free energy equation as follows:
(III.12)
Throughout this equation and referring to (Harvie et al. 1984; Pitzer 1979, 1991), the subscripts M and b refer to cations, similary X and a refer to anions. The values ZM, ZX, mb and ma are the charge and molal concentrations of cations and anions respectively. , refer to the summation of the proprieties over all anions and all cations respectively. , denote the sum over all possible pairs of variation of all cations and anions respectively (Boris 2001)
The salt dissolution-precipitation progress rate vs is expressed by a Fick’s law as follows:
(III.13)
The kinetic constant ks[mol.T-1] depends on the thermodynamic conditions, the species concentrations and on the reactant surface extent.
A kinetic approach is used in the present model, at each time step after calculating the activities of all species in the solution, a test on the equilibrium state of each salt, existing in the system or that could be generated, is performed. The precipitation or the dissolution of this salt is quantified according to the progress rate vs. Regarding the high nonlinearity of the mathematical formulation due to the ion activity calculation according to the Pitzer model, a successive iteration is carried out between the dissolved or precipitated salts quantities and the solution species concentrations that determine their activities. Obviously, the Gibbs phase rule (W. Stumm and Morgan 1970) has to be verified (the phase rule is an expression that relates the number of independent components (C), the number of phases (P), and the number of degrees of freedom (F) in an equilibrium system), as well as the electrical neutrality of the solution. Subsequently, for the reached composition, a precise value of the solution density is calculated according to the mixture excess free energy formulation of Pitzer (Monnin 1989). Ten species are retained at this stage of the model development: ,,,, ,, , , and , and a list of nineteen evaporitic salts that can be generated in the brine evolving in a neutral saline way (Bouhlila and Laabidi 2008b). The database developed in this work contains for each species, the chemical formula, the molar mass, the partial molar volume, the charge number and the Pitzer coefficients of the species for the chemical activities and the solution density calculations. It contains also for each salts the solubility product, the precipitation and the dissolution kinetic constants, the stœchiometric coefficients, the molar mass and the molar volume.
Obviously, each particular case to be processed is described by its specific data such as the initial composition of the solution, the reason of its equilibrium displacement: by evaporation, by mixing of two waters or by reaction with new mineral phases, and the total simulation time
Due to the high nonlinearity of the activity expressions according to the Pitzer model, the following approach was adopted in the dissolution-precipitation reactions:
At each time step t and after calculating the activity of all aqueous species, the ion activity product of each salt IAPs is calculated and ranked in ascending order;
Dissolution: the salt with the lowest saturation index is dissolved if it is available in enough quantity in the solid matrix according to the kinetic law. The same protocol will be restarted for another salt with a higher saturation index, until no further dissolution is possible.
Precipitation: after updating the activities of all the species (after the dissolution reactions), the saturation index is ranked in ascending order, after that we gradually precipitate the salts with a saturation index greater than unity starting by the highest one and according to the kinetic law.
III.III. Hydrogeochemical model: Coupled model
III.III.1. Governing equations
GEODENS describes transient convective, diffusive-dispersive and reactive transfer of specie (=1,) using a sequential iteration approach with the following equation:
. (III.14)
[MM-3] is the concentration of solute, and is the total number of species, and represents the transfer operator applied to including the convection and the diffusion-dispersion phenomena and is expressed as follows:
. (III.15)
The term represents the geochemical mass flux of species. If Nr is the number of reactions related to species , is calculated as the sum of the reaction contributions into this mass flux (Bouhlila 1999a), and is written as follows:
(III.16)
The geochemical flux of each species at Δt is calculated as the mass change of this species during dissolution-precipitation reactions.
III.III.2. Existing permeability–porosity relationships and implementation into reactive transport models
The physicochemical processes are written at the macroscopic scale in GEODENS. The pore geometry is continuously modified by the dissolution-precipitation reactions of salts. The porosity at any node of the mesh and at every time step Δt is updated as follows:
(III.17)
where Vjs [L3M-1] is the molar volume of the salt j, rjs[M-1L-3T-1] is the rate of the mineral precipitation (-) or dissolution (+), and Ns is the total number of salts.
(Le Gallo et al. 1998) presented a brief review of porosity-permeability relationship models based mainly on reservoir engineering literature. Model applications focused especially on the influence of reaction around well bores. For the modeling of the relation between porosity and permeability, capillary tubes or plane cracks models have been used and can give relatively simple relationships between these two parameters such as the Carman-Kozeny model (De Marsily 1981b). In the following we quote the popular porosity-permeability relationship available in the literature:
Brinkman model: the Brinkman model considers flow around a sphere, and links the drag on the fluid to permeability by assuming the pore walls are barriers to the fluid flow and is written as follows:
(III.18)
Carman-Kozeny model: Carman-Kozeny model is defined as follows:
(III.19)
where is the tortuosity of porous media and S is the specific surface area.
Fair-Hatch model: A modified Carman-Kozeny model is provided by Fair-Hatch model (Le Gallo et al. 1998):
(III.20)
where is the mineral volume fraction, is the grain radius, and is the number of reactive minerals.
(Sanford and Konikow 1989a) proposed an empirical model for the sedimentary rocks to represent the permeability variations due to the dissolution of calcite in a coastal aquifer as follows:
(III.21)
According to these authors,a is about 20 and b is deduced from the initial known values of k and.
Model for fractured media: A simple model for fractured media was presented by (C. I. Steefel and Lasaga 1994) based on case of a set of parallel fractures with smooth walls and is written as follows:
(III.22)
where is the fracture porosity, is the fracture aperture, and d is the fracture spacing.
IV.I. Integration of coupled equations
In the integration of the flow and transport equations, the flow domain, Ω is discretized by using triangular finite elements, Ωe with six nodes, n1, n2, n3, n4, n5 and n6 as shown in Fig. IV.1.
If we consider that ν=(ν1, ν2), the output unit normal to the boundary line of Ω, the boundary conditions can be expressed as follows:
For the flow: In Г1 the pressure (p=p1) is imposed and in Г2 a flux (q=q1) is imposed.
For the transport: In Г3 the concentration (C=C1) is imposed and in Г4 a convective or dispersive-convective flux () is imposed.
And the following relationship is verified:
(IV.1)
Figure IV.1 Flow domain and subflow domain (triangular finite elements)
The Galerkine approximation method is used to discretise the main field variables, p(x,y,t) and C(x,y,t). The error (residual) generated by the introduction of this nodal approximation in the governing equations needs to be minimized over the entire flow domain. The weighted residual method is based on the elimination of this residual. To achieve this goal the residual is weighted, based on an appropriate function and is integrated over the flow domain Ω as follows:
(IV.2)
Figure IV.2 Transformation reference element and real element
The nodal approximation of the main field variables are as follows:
(IV.3)
are the interpolation functions in node I, I=1,n.
are calculated as follow (Dhatt and Touzot 1981).
1 if I=J
= (IV.4)
0 if no
The integration by parts of the flow and transport equations leads to a lower order integral form. By this method a matrix form of the problem can be formulated (which is a differential equation problem).
By the application of Galerkine method and by using integration by parts of the fluid mass conservation equation,
(IV.5)
the following low order integral is obtained:
(IV.6)
the matrix form of this integral is expressed as follows:
(IV.7)
Integration by parts of the nonconservative form of the transport equation is carried out as follows:
(IV.8)
integration by parts of the diffusive term only gives the convective form of the transport equation CFTE (H. J. Diersch 1988; Giorgio Galeati and Gambolati 1989).
(IV.9)
the matrix form of this integral is expressed as follows:
(IV.10)
where is the diffusive flux exchange at the node I.
By application of the Galerkine method to the conservative form of the transport equation and through integration by parts of the diffusive and advective terms we obtain the diffusive form of the transport equation, DFTE:
(IV.11)
the matrix form is expressed as follows:
(IV.12)
where is the total flux (convective, diffusive and dispersive) exchange at the node I.
The convective form (CFTE) of the transport equation leads to a more stable numerical scheme, especially when convection is important (Bouhlila 1996; H. J. Diersch 1988; Giorgio Galeati and Gambolati 1989).
The Darcy velocity components are used in the transport equation, for the convective term, and for the dispersion tensor calculation. The approach used to evaluate the Darcy velocities in this work is called the discontinuity velocity approach. It consists of calculating the velocity in each element of the mesh from the pressure and concentration values in the nodes by a simple application of Darcy law. Obviously, these velocities are discontinuous between elements which can generate numerical instability especially in the cases where advection is significant (H. J. Diersch 1988). In fact, and as the calculation of velocity incorporates the difference between and,it is necessary that these two terms have a similar order. To avoid erroneous vertical velocity components, Voss (1984) adopts a consistent balancing of the gravity term by using the derivatives of the nodal interpolation functions on quadrilateral elements of four nodes. In this work a different method is adopted which leads to the same result for the variable approximation to the Darcy velocity. It consists of calculating the gravity term for a triangular element of six nodes, by using the concentration values of the three summit nodes and applying linear nodal interpolation functions.
After the finite element discretization of the variable density flow and transport, further inputs to the model code which are related to the nonlinear and to the nonstationary algebraic systems are needed. The time discretization scheme is quasi-implicit. The first derivate of time is approximated at time (t +αΔt), where using a finite difference expression:
(IV.13)
where Xt and Xt+Δt , respectively, the values of the variable X at time t and t+Δt respectively. The value of X at (t +αΔt) is calculated as follows:
(IV.14)
An iterative calculation is used in order to evaluate all the matrices of these nonlinear systems at t+Δt.
During the first iteration a linear extrapolation of p and C is used as a first approximation of these variables. These extrapolation values {p}ext and {C}ext are evaluated at t+Δtn+1 as follows:
(IV.15)
These values are used to start the first iteration in solving the nonlinear matrix system. On further iterations the system is updated by using the values and of the previous ones. These values are used also in the calculation of the Darcy velocities in every element of the mesh. The iterative calculation at t+Δtn+1 is finished when the following convergence criteria for p and C are achieved:
(IV.16)
whereandare the difference between two successive Iterations of p and C respectively. The update of the variable vector (p and C) is obtained via a relaxation coefficient ω as follows:
(IV.17)
The hydrogeochemical model is the coupling between the density dependent flow and the transport model extended to ten concentrations with the geochemical model where a sequential iterations method is retained.
By the application of the Galerkine method only on the geochemical flux, we obtain:
(IV.18)
After incorporating the geochemical flux term in the CFTE form we obtain for every aqueous species i the following matrix system:
(IV.19)
IV.II. Calculation strategy and resolution
As shown in flow chart presented in Figure IV.3, the resolution is mainly implemented by using two steps: the input reading and the transient simulation.
Input reading: in the hydrogeochemical modelling the inputs are classified into hydrodynamic and geochemical data.
-The hydrodynamic data: The finite element mesh containing the elements and nodes list, element proprieties (permeability, porosity, dispersivity, and diffusion coefficient), boundary conditions (boundary conditions values in pressure and concentrations of all the species), initial conditions (initial conditions in pressure, in concentration of the species, and in the quantity of salts in each elements of the mesh), time steps and simulation time.
– The geochemical data: as described before (Chapter III) the required data is mainly the data related to the ions, salts and the Pitzer coefficients.
Transient simulation and geochemical flux calculation: the method used is explicit using the concentrations values of the previous iteration.
Using the concentrations Ci(t-t), the chemical activities of aqueous species i are calculated, then the solubility products of the selected salt is tested. The dissolution-precipitation is then quantified by using a kinetic or a chemical equilibrium approach.
In the case of chemical equilibrium approach, the concentrations C of this equilibrium are calculated and are related to the geochemical flux of the specie i at the node I as follows:
(IV.20)
In fact, due to the nonlinearities in terms of activities, the calculation of equilibrium concentrations C is not straightforward especially when there are several species in solution and many reactions are considered. The dissolved or precipitated quantity of salts is calculated using a successive iteration method. The equilibrium is considered achieved when all are close to 1.
In the case of the kinetic approach, the geochemical fluxes of the specie i at the node I is obtained as follows:
(IV.21)
where Nr is the number of the chemical reaction.
Figure IV.3 Flow chart of the numerical procedure in GEODENS
The choice between these two approaches is determined by comparing the time step t and the kinetic coefficient Kcr of each reaction r as follows:
If the reaction r reached the equilibrium before the end of the time step Δt, so the equilibrium approach is considered.
if the reaction r is partially progressed during the time step Δt, so the kinetic approach is considered.
CHAPTER 5. MODEL VERIFICATION
Analytical solutions are not available for complex multicomponent reactive transport problems. Previous laboratory study results have been employed to overcome this shortfall. The verification of the present model consists firstly to check the density dependent flow module, secondly to check the geochemical module and finally, to verify the results of the coupled model. Regarding the density dependent flow model, the first benchmark used is the Henry problem, (H. R. Henry 1964a) solved the governing equations of the density-dependent flow in a non-dimensional form and presented a semi-analytical solution for the steady state distribution of salt concentration and stream functions. The second benchmark used is the convection cell problem which describes the free-convection induced by buoyancy effects. The geochemical component of the GEODENS code was successfully verified using results on the effect of the saltwater fraction on the chemical processes (dissolution-precipitation reactions) in a carbonate environment during the mixing between saltwater and freshwater, as described further below. The validation dataset was presented by (Rezaei et al. 2005). And finally, the coupled model was validated by using the laboratory presented by (Singurindy et al. 2004a), a simulations of the porosity evolution in a 2D flow cell due to calcite precipitation-dissolution reactions was carried out and compared to the laboratory experiments results.
V.I. Benchmarking Density dependent flow model
V.I.1. Rayleigh problem
The density dependent flow and coupled transport model with the different numerical schemes developed are first verified against the cell density driven problem. Although it is a simplified situation, the related results can shed light on other more complex and more real configurations. This cell problem can also be used as a test case for density-dependent groundwater flow and solute transport simulators (Weatherill et al. 2004).
The process in the convection cell is governed by the Rayleigh number noted Ra. Ra is a measure of the convective stability and is represented as the ratio of buoyancy forces over the resistive forces. Buoyancy forces give rise to a free-convection flow that in turn increases the dispersion and breaks down density contrasts. Ra is calculated as follows:
(V.1)
Where: is the difference between the density in the upper and in the lower boundary, μ is the fluid viscosity, k is the permeability, is the porosity and d0 is the diffusion coefficient. As presented in (Elder 1967; Scheidegger 1974) three conditions were presented:
– if Ra4π2: there is no convection flow
– if 4π2Ra240 to 390 : stable convection flow
– if Ra240 to 390 : unstable convection flow.
Besides, multiple solutions for this flow density driven problem can exist and only one is obtained during a given time-evolution simulation (H. J. G. Diersch and Kolditz 2002; D. Henry et al. 2012). As an element for comparison between the solutions, and also to quantify the contribution of convection to the transport, we introduce the Sherwood number noted Sh which is the ratio of the total flux of the salt (sum of the convection and the diffusive-dispersion flux) by the diffusive flux. She is defined as follows (H. J. Diersch 1988; Elder 1967):
(V.2)
The geometry for the convection cell is shown in Fig. V.1 and the related physical parameters are summarized in Table V.1. The boundary conditions imposed in this problem are shown in Fig. V.1 where a no-flow boundary is applied to all the domain, a concentration C0 is applied to the lower boundary and a concentration C1>C0 is applied to the upper boundary. Two finite elements meshes, of 441 nodes and 121 nodes respectively, are used in these comparative simulations.
Table V.1 Model parameters used in the convection cell simulations
Figure V.1 Convection cell domain and boundary conditions.
Numerical schemes comparative study :
We summarize in this section the results for the saline Rayleigh cell simulated with the two forms of the transport equation, DFTE and CFTE, discretized with the described Galerkine finite elements method, GFEM.
For all scheme, the presented results correspond to a quasi-stationary regime reached with transient simulations starting from a hydrostatic state where the entire cell is at the concentration C0, the concentration C1 is applied to the cell top boundary, starting at time 0. The results for the different simulations compared to the previous studies are shown in Table V.2. We compare the Sherwood number calculated for the quasi-stationary regime reached in the two numerical schemes with those calculated by (Caltagirone 1975; H. J. Diersch 1988). It appears that the GFEM-CFTE scheme, with the refined mesh, provides the most satisfactory results.
The concentration distribution and the velocity field for different Ra number are shown in Fig. V.2.
Table V.2 Sherwood number calculated by different numerical schemes and compared with those of (Caltagirone 1975; H. J. Diersch 1988)
(-) the numerical scheme doesn’t converge * Unstable convection flow
Ra: Rayleigh number Sh : Sherwood number
GFEM: the Galerkine finite element method C(A): the calculated concentration at A(x=0.5,z=0.8)
DFTE: the diffusive form of the transport equation CFTE: the convective form of the transport equation
Figure V.2 Concentration distribution and velocity field for different Rayleigh number: Ra=100 (a), Ra=250 (b), and Ra=500 (c).
V.I.2. Henry problem
Henry problem is widely used as a first test of accuracy during the development of a variable-density code. Henry problem has been studied extensively (Ackerer et al. 1999; Benson et al. 1998; Boufadel et al. 1999; Croucher and O'Sullivan 1995; H. R. Henry 1964a; Oldenburg and Pruess 1995; Putti and Paniconi 1995; C. I. Voss and Souza 1987; Weatherill et al. 2004). It concerns seawater intrusion into a coastal aquifer. A rectangular domain represents a vertical cross-section: its left-hand side has a constant inflow of freshwater, while on the right there is a seawater boundary (Fig. V.3). The top and bottom boundaries are impermeable. Concentration is scaled so that c = 1 corresponds to seawater concentration. The parameters given in Table V.3 are those which match Henry's original paper,
Figure V.3 The Geometry and boundary conditions of the Henry problem.
Table V.3 Parameters used in Henry’s problem.
The problem requires only a steady-state solution (C. I. Voss and Souza 1987). It involves three out of the four physical processes which may be involved in variable- density flow and transport: regional flow, density-induced convection and diffusion but not dispersion.
Figure V.4 Henry problem: Concentration contours and velocity field after 100 minutes.
Fig.V.4 displays the concentration contours and velocity field after 100 minutes as (C. I. Voss and Souza 1987) run SUTRA for 100 minutes with time steps of 1 minute, and assume that the simulation has reached steady-state. However, for codes using the same boundary conditions, there is only minor disagreement between this work and other codes such as SUTRA (C. I. Voss and Souza 1987) and CODESA3D (Gambolati et al. 1999) as shown in Fig. V.5.
Figure V.5 Henry problem: Concentration distribution at domain bottom compared to the result of CODESA3D and SUTRA.
V.II. Benchmarking geochemical model
The effect of the saltwater fraction on the dissolution-precipitation reactions in the carbonate environment during the mixing between a seawater and freshwater were discussed and presented by many authors. (Rezaei et al. 2005) described the mixing process by a 1D-model representing a mixing line, with end members concentrations prescribed at either side respectively.
In the present work, a simple geometry has been used to simulate the geochemical processes in the mixing of saltwater and freshwater equilibrated with calcite (solution 1 and solution 2 respectively, Table V.4).
Table V.4 Chemical composition of boundary solution used in the simulation.
Units for (total) concentrations are g Kg-1
Mixed waters are under-saturated with respect to calcite throughout the domain during the transport of the chemical species by diffusion. Fig. V.6a shows the cumulative quantity of dissolved calcite during the mixing process in the calcite solid matrix and Fig. V.6b displays the calcite saturation index. Fig. V.6 shows that maximum dissolution occurs much closer to the freshwater side at salinities below 20%. A similar result is also obtained analytically by (De Simoni et al. 2005a). This has also been found by numerical modeling (Romanov and Dreybrodt 2006; Sanford and Konikow 1989a).
Figure V.6 Mixing of freshwater and saltwater equilibrated with calcite (solutions 2 and 1, respectively, Table V.3): (a) calcite saturation index; and (b) cumulative quantity of dissolved calcite (expressed in mol/L) of mixture in a carbonate reactive medium.
We can conclude that there is a good agreement between the equilibrium method used by (Rezaei et al. 2005) and the method used in this work. In summary the chemical equilibrium method used in GEODENS is satisfactory.
V.III. Benchmarking the coupled model
(Singurindy et al. 2004a) conducted a laboratory experiment of the carbonate dissolution-precipitation process in a 2D cell during mixing between freshwater and saltwater. These laboratory experiments focused on the porosity changes due to the precipitation or the dissolution of the carbonate rocks during the mixing process.
Referring to (L. N. Plummer et al. 1979), the two most important factors generating the precipitation or the dissolution of such solids are the amount of seawater in the mixture and the CO2 partial pressure. If the fraction of seawater in the mixture zone is low or the CO2 partial pressure is high, then these two conditions can generate the dissolution of calcite in the groundwater. A real example of carbonate dissolution is found in coastal aquifers such as the Yucatan Peninsula, Mexico (Back et al. 1979a). If the fraction of seawater in the mixture zone is high and PCO2 is low, then the precipitation of calcite in the groundwater can be generated. These factors could explain the absence of dissolution in the coastal aquifer of Mallorca, Spain (Price and Herman 1991a).
An analytical solution to calculate the porosity change in the mixing zone was suggested by (Phillips 1991). The feature of this result is the decoupling of the transport and of the geochemical reaction, because if dissolution of the mineral is sufficiently fast, then the solution can be regarded to be in equilibrium with respect to the dissolving mineral. These equations can be decoupled and the computational effort is highly reduced (Romanov and Dreybrodt 2006).The change of the porosity is caused by the mixing of the freshwater and saltwater that is in equilibrium with respect to calcite and under a constant CO2 partial pressure. The rate of the porosity change is proportional to the macroscopic dispersion coefficient, to the second derivative of the equilibrium concentration of the mixture, and to the square of the salinity gradient (Phillips 1991).
Figure V.7 Geometry and boundaries conditions of the modeled domain used in the Laboratory experiments.
In order to validate our approach, we used GEODENS to simulate the flow, transport and geochemical reactions in the geometry used by (Singurindy et al. 2004a) in their laboratory experiment. The basic domain is shown in Fig. V.7. The domain is 16 cm×16 cm×1 cm; two reservoirs supply the domains, one containing a freshwater solution and the other containing an NaCl solution. The two waters were saturated with calcite and under atmospheric pressure PCO2=1atm for the dissolution experiment and at PCO2=10-3.5 atm for the precipitation experiment.
Table V.5 Chemical composition of boundary solution used in the simulations.
Units for (total) concentrations are g Kg-1
A grid of 441 elements was used to simulate this laboratory experiment.
V.III.1. Dissolution experiment
The first reservoir (pure water solution) supplies the domain through the upper boundary under a fixed head boundary (h=1cm). The second reservoir supplies the domain through the right hand side by an NaCl solution of 40 g/l under the same fixed head value. In Table V.5, we present the chemical composition of the solutions at the two fixed heads boundaries of the domain used in the calculations. The two remaining boundaries present the outlets of the domain under fixed head boundary conditions equal to zero (h=0m). A dispersivity of 6.5 x 10-4 m and a porosity of 0.33 are used.
Fig. V.8 shows the calculated concentration distribution of the chloride (Cl-) in the modeled domain after the mixing of the two solutions for dissolution experiment.
Figure V.8. Cl- isolines illustrate the calculated salinity distribution for dissolution experiment
Fig. V.9 shows the spatial distribution of the porosity change for the dissolution experiment. We have found the same structure and location of the porosity change as in the simulation of (Romanov and Dreybrodt 2006) and the laboratory experiment of (Singurindy et al. 2004a). The laboratory experiment shows a porosity change along the domain diagonal (the mixing zone) as in the present work.
Figure V.9 Spatial distribution of porosity change for dissolution experiment after 6 hours.
The change of porosity has been determined as 0.240×10-4 year-1 for dissolution in the laboratory experiment of (Singurindy et al. 2004a).In our work, an average of 0.238×10-4 year-1, this value is reached after many simulations using different Kc values as presented in Table V.6. The results show that the simulation with a kinetic coefficient value equal to 1.15×10-10 moles/s is in good agreement with the experiment. In Table V.6, we present the dissolved quantity of calcite and the change of porosity evaluated analytically, experimentally and numerically (this study).
Table V.6 Theoretical, experimental and numerical (this study) estimates of dissolution rates and porosity change.
The results show that the approach presented in this work is in good agreement with the experiment and the previous numerical studies.
V.III.2. Precipitation experiment
The first reservoir (pure water solution) supplies the domain through the upper boundary under a fixed head boundary (h=1cm). The second reservoir supplies the domain through the right hand side by an NaCl solution of 80 g/l under the same fixed head value. In Table V.5, we present the chemical composition of the solutions at the two fixed heads boundaries of the domain used in the calculations. The two remaining boundaries present the outlets of the domain under fixed head boundary conditions equal to zero (h=0m). A dispersivity of 6.5 x 10-4 m and a porosity of 0.33 are used.
Fig. V.10 shows the calculated concentration distribution of the chloride (Cl-) in the modeled domain after the mixing of the two solutions for precipitation experiment.
Figure V.10 Cl- isolines illustrate the calculated salinity distribution for precipitation experiment
Fig. V.11 shows the spatial distribution of the porosity change for the precipitation experiment. We have found the same structure and location of the porosity change as in the simulation of (Romanov and Dreybrodt 2006) and the laboratory experiment of (Singurindy et al. 2004a). The laboratory experiment shows a porosity change along the domain diagonal (the mixing zone) as in the present work.
Figure V.11 Spatial distribution of porosity change for precipitation experiment after 6 hours.
The change of porosity has been determined as -0.05 year-1 for the laboratory experiment of (Singurindy et al. 2004a).In our work, an average of -0.05 year-1, this value is reached after many simulations using different Kc values as presented in Table V.7. The results show that the simulation with a kinetic coefficient value equal to 2.15×10-12 moles/s is in good agreement with the experiment. In Table V.7, we present the precipitated quantity of calcite and the change of porosity evaluated analytically, experimentally and numerically (this study).
Table V.7 Theoretical, experimental and numerical (this study) estimates of precipitation rates and porosity change.
The results show that the approach presented in this work is in good agreement with the experiment and the previous numerical studies.
VI.I Model Application
In this part of this chapter, we check the capability of GEODENS code (Bouhlila 1999b; Bouhlila and Laabidi 2008a; Laabidi and Bouhlila 2015) in calcite dissolution and precipitation reactions and their effect on the flow and transport problem. The model tested and validated is used to simulate two applications. The first application treated the clogging effect on the flow and transport problem during calcite precipitation reaction, and in order to check the capability of GEODENS code (Bouhlila 1999b; Bouhlila and Laabidi 2008a; Laabidi and Bouhlila 2015) in precipitation and deposition reaction, the same experiment conducted by (Katz et al. 2011) was carried out in this work. The second application treated the effect of the calcite dissolution on the penetration length of the seawater intrusion, and in order to provide general results we simulate these processes using Henry problem geometry with the same hydrodynamic parameters and boundary conditions.
VI.I.1. Precipitation reaction simulation
Permeability reduction is an important process that needs to be considered especially during groundwater injection. The concern associated with porosity and permeability reduction due to calcite precipitation is identified in many environmental and industrial applications including CO2 geological injection and in situ groundwater remediation. In the following the impact of calcite precipitation on these two techniques were discussed briefly:
During CO2 geological injection, carbonates are the first minerals to react. Along the pathway of the injected CO2 during geological sequestration, calcite is formulated by precipitation reaction. The calcite deposition can easily decrease the porosity and permeability of the solid matrix that can reduce reservoir capacity and decrease the injection efficiency and the available storage volume (Zhang et al. 2010).
As reported by Fujita et al. (2008) during groundwater remediation subsurface microorganism can generate calcite precipitation by increasing pH of surrounding fluids. This phenomenon may lead to a decrease of the aquifer porosity and permeability or to a reduction of the reactive surface of the treatment system. Calcite deposition by precipitation reaction along the plume margin during in situ remediation can reduce the cleanup efficiency (Zhang et al. 2010).
For a given mineral, if the ion product is greater than the solubility product, then the solution is supersaturated with respect to this mineral and precipitation can occurs (Werner Stumm and Morgan 1981). Supersaturation may be occurred by three main mechanisms: (i) changes of fluid pressure or/and temperature during flow through porous media, (ii) dissolution of a given mineral that results a supersaturation with respect to another one (e.g. dissolution of calcite can induces precipitation of gypsm (Singurindy and Berkowitz 2003b)), and (iii) mixing induced precipitation. Referring to (L. N. Plummer 1975b) the two most important factors generating the precipitation or the dissolution of such solids are the amount of the seawater in the mixture and the partial pressure of carbon dioxide pCO2. If the fraction of the seawater in the mixture zone is high and pCO2 is low then, the precipitation of the calcite in the groundwater can be generated. These factors could explain the absence of dissolution in the coastal aquifer of Mallorca, Spain (Price and Herman 1991b). Mineral deposition in geological formation and decrease of porosity due to the infilling of fractured and porous media during precipitation reactions are of great interest in many environmental and industrial applications such as effluent disposal (Eriksson and Destouni 1997) and oil recovery extraction (Rege and Fogler 1989). The lack of data on precipitation kinetic in subsurface is the main cause of the lack of understanding of calcite filling in porous media. Modeling of expected time and amount of calcite needed to produce mineral filling has been a largely unresolved problem.
In order to check the capability of GEODENS code (Bouhlila 1999; Bouhlila and Laabidi 2008; Laabidi and Bouhlila 2015) in precipitation and deposition reaction, the same experiment conducted by Katz et al. (2011) was carried out in this work. The geometry and boundary condition was illustrated in Fig. 4. The flow cell dimension used is 25 cmx10 cmx0.8 cm in the x-, the y- and the z-directions respectively saturated with an initial aqueous solution of sodium chloride. A fixed head boundary is applied along the inlet face and the outlet face is open to flow at head zero, such that an average of water velocity of 0.97 m/day inside the flow cell is created. The inlet face is modeled in such a way to create two separate inlets, side by side, each with a length of 5 cm. Initial value of permeability was estimated by Katz et al. (2011) using kozeny-Carman model and was set to a value k0=1.12×10-9 m2 and the initial porosity to a value Φ0=0.375. for chemical transport, it was assumed that all chemical species are characterized by the same molecular diffusion D0= 10-9 m2/s. the transverse dispersion has an important role because it governs mixing process in the case of parallel flow and is reached after a best fit of Ca2+ concentration.
Conservative experiment simulation
The aim of this simulation is to calibrate the transverse dispersivity, we assume in the case of a parallel flow that longitudinal dispersivity hasn’t any effect on the flow and transport problem. Fig.VI.1 depicts the flow cell geometry and the boundary conditions. A grid of 441 elements is used in this simulation and the hydrodynamic proprieties assigned are described above. Two solutions of sodium chloride NaCl and calcium chloride CaCl2 are injected in parallel.
Figure.VI.1. Flow cell geometry and boundary conditions for the conservative transport simulation.
In order to fit the numerical simulation results to the laboratory experiment results (Katz et al. 2011), we proceed to measure the Ca2+ concentrations in some sampling ports. Legend and locations of these sampling ports are given in Table VI.1 and illustrated in Fig. VI.1.
Table VI.1. Coordinates of sampling ports used to measure the Ca2+ concentration (Precipitation reaction simulation).
αT is changed several time in order to get a best fit between numerical and laboratory experiment results in sampling ports. Comparison against measured time evolution of Ca2+ concentration resulted in an overall best fit value αT =0.2 mm. Fig. VI.2 shows comparison between model results (solid line) and laboratory experiment results (symbol) conducted by Katz et al. (2011) at the selected sampling ports. Our results seem to be in good agreement and the model capture the spatial and temporal distribution of the Ca2+ concentration.
Figure.VI.2. Calibration of dispersivity. Total aqueous calcium concentrations in conservative transport simulation (for model prediction (solid line) and experiment results of Katz et al. (2011) (symbol)), at locations of the sampling ports.
From Fig. VI.2a,c,e, and g we can depict that the model predicts reasonably accurate concentration in port located at calcium side (CaCl2 boundary) however Fig.VI.2 b,d,f, and h shows that the model capture the Ca2+ distribution as in laboratory experiment results but with a higher magnitude this is caused basically by the location of these ports, to the non calcium side (NaCl boundary) and due also to the characteristic of the calibration exercise, where a change in the despersivity value would improve the fit in a ports and degrade that in the other.
Reactive experiment simulation
The aim of this simulation is to validate our coupled model regarding the laboratory experiment results conducted by (Katz et al. 2011) and to check the capability of GEODENS code (Bouhlila 1999b; Bouhlila and Laabidi 2008a; Laabidi and Bouhlila 2015) on precipitation and deposition reaction during mixing. In this regard, the same geometry and hydrodynamic proprieties as in the conservative experiment simulation were used. The inlet boundary is exposed to two solutions Na2CO3 and CaCl2 injected in parallel as shown in Fig.VI.3.
Figure.VI.3. Flow cell geometry and boundary conditions for the reactive transport simulation.
The fitted value αT =0.2mm is used in this simulation. Calcium carbonate precipitated along the mixing zone in the cell centerline between the two solutions Na2CO3 and CaCl2. In this simulation the precipitation reaction is modeled as kinetic reaction and a value of kinetic coefficient kc of 7.5×10-10 moles/s is used.
As shown in Fig.VI.4 even the first deposition is created near the intlet boundary, fluids are able to flow around closed regions and into more permeable area where mixing can occur in the other parts of the cell centerline. This process is repeated until all the mixing region is filled by the precipitated mineral.
Fig.VI.4a and Fig.VI.4b display respectively the precipitated calcium carbonate amount and the porosity change within the flow cell at 60h. Comparing the results presented by (Katz et al. 2011) where they found that the large reduction of porosity is located near to the inlets boundary, our model calculates the same magnitude of porosity change as presented by (Katz et al. 2011) but located along the centerline cell. A high amount of precipitated calcium carbonate is located on the mixing zone which is about 1×104 mol/m3 and a reduction of porosity reaching a value of 8×10-2.
Fig.VI.5 shows the temporal distribution of Ca2+ concentration during simulation. Generally during mixing inducing precipitation of calcium carbonate, we expect a decrease of Ca2+ concentration. This can be seen clearly by comparing the measured concentration at sampling ports in the reactive simulation to the corresponding values detected for the conservative one.
Fig.VI. 5b,d,f, and h show the Ca2+ concentration at the ports A3, A4, C3 and D3 which are located in the Na2CO3 side. As discussed in (Katz et al. 2011) the sharp decrease of Ca2+ concentration at sampling ports may be associated to the effect of the clogging process. During the first simulation hours, clogging process has a direct effect on the flow and transport problem, such as a decrease of Ca2+ (Tartakovsky et al. 2008) concentration, that can be identified in these sampling ports. Concentration continued a decrease until a minimum value of 5×10-5 mol/kg-water found in port D3. Referring to (Tartakovsky et al. 2008) if clogging process separates the two solutions, then they may become undersaturated with respect to calcite and dissolution can occur. In this case it is possible that an increase of Ca2+ concentration can be identified at later stage of the simulation. On the other side of the flow cell, in the CaCl2 boundary we can conclude that our results are in good agreement with the laboratory experiment results. The values recorded in these ports are significantly lower than those in conservative scenario. This is associated as described before especially with the fact that Ca2+ is subtracted from the solution during mixing inducing precipitation of calcium carbonate.
Comparing to the results presented by (Katz et al. 2011), we can conclude that GEODENS code (Bouhlila 1999b; Bouhlila and Laabidi 2008a; Laabidi and Bouhlila 2015) appears to capture the majority of the processes that can be identified in the laboratory experiment but with different magnitude in Ca2+ concentration, in the amount of precipitated calcite and in the porosity change, and this can be explained by the precipitation rate used in this simulation.
Figure.VI.4. Simulation results of reactive transport after 60 h. (a) Calcite accumulation; (b) spatial distribution of porosity.
Figure.VI. 5. Total aqueous calcium concentrations in reactive transport simulation (for model prediction (solid line) and experiment results of Katz et al. (2011) (symbol)), at locations of the sampling ports.
VI.I.2. Dissolution reaction simulation: Reactive Henry problem
As reported by Cabral et al. (1992) and Bear and Kapuler (1981) , during seawater intrusion the flow is assumed to be perpendicular to the short line, so that two dimensional vertical model is satisfactory to model such problem. In this context, Bear (1979) has given some examples of real aquifers in which field data have shown that saltwater concentration is distributed in a transition zone sharply separating freshwater and saltwater. His main assumption is to neglect the transition zone compared to the aquifer dimensions. In this study, the same geometry and parameters as in Henry problem are used (Henry 1964). The Henry saltwater intrusion problem concerns a vertical slice through an isotropic, homogeneous, confined aquifer. For the flow boundary, a freshwater constant flux is applied from the inland boundary and a hydrostatic pressure is distributed on the sea boundary. As shown in Fig. 1 it is not easy to represent the transport boundary in the seaside. In this work a mixed boundary condition in the seaside is introduced in order to avoid difficulties in the finite element solution near the outflow region of the aquifer (Segol et al. 1975). Segol et al. (1975) replaced the original homogeneous Dirichlet boundary by a mixed Dirichlet-Neumann boundary. The upper 20 cm was a zero concentration flux exit boundary, and the lower 80 cm remained a constant concentration Dirichlet boundary.
Figure.VI.6. Henry problem: domain and boundary conditions, Henry (1964).
Henry (1964) solved the governing equations in a nondimensional form and presented a semi-analytic solution for the steady state distribution of salt concentration and stream function. Henry problem has been widely used as a test case (benchmark) of density dependent groundwater flow models (Abarca et al. 2007; Simpson and Clement 2003; Voss and Souza 1987). The boundary conditions of Henry problem are defined in Fig.VI.6. No flow occurs across the top and bottom boundaries. A freshwater flux is injected along the left vertical boundary of the domain. The prescribed pressure along the right vertical boundary is set at hydrostatic seawater pressure. A grid of 231 elements is used to simulate this problem. The soil parameters used in Henry problem are summarized in Table VI.2. In the original Henry problem, a pseudo-permanent regime is reached after 100 minutes (Bouhlila 1999; Voss 1984). If the solid matrix is calcareous, the dissolution processes are perpetual and no permanent regime is reached.
Table VI.2. Parameters used in Henry’s problem.
Recently, Henry reactive problem was presented by Laabidi and Bouhlila (2015) where the authors used an alternative modeling approach to simulate the effect of calcite dissolution on the porosity and permeability. The authors showed through their simulations that porosity change due to calcite dissolution has a direct effect on the penetration length of seawater and the flux from the seaside.
The results presented by Abarca et al. (2007) and by Nick et al. (2013) showed that toe penetration of saltwater decreases with increasing transversal dispersivity, while width of the mixing zone and saltwater flux become larger.
In this simulation, the effect of calcite dissolution reaction on the behavior of saltwater intrusion during mixing is tested. During dissolution, porosity and permeability change is identified and the consequence of this latter on the flow and transport is calculated. Several simulations are carried out in order to get a better understanding of the effect of calcite dissolution on seawater intrusion using different kinetic coefficient values ranging from kc=10-11 moles/s to kc=9×10-10 moles/s. Fig. VI.7A to Fig.VI.12A show the porosity change for different kc values due to calcite dissolution during the mixing of the saltwater and freshwater after a simulation time of 25, 50, and 100 years. The quantity of the dissolved calcite is extremely related to the rate of calcite dissolution and to the order of the mixing (fraction of the saltwater). Maximum values are located in the mixing zone (transition zone) between seawater and freshwater. Fig.VI.7B to Fig.VI.12B show the permeability change for different kc values due to porosity change. As the porosity change, maximum values are located in the transition zone between freshwater and seawater. This new hydrodynamic distribution (porosity and permeability) controls the variation of the velocity field and the shape of the concentration repartition as shown in Fig.VI.7C to Fig.VI.12C which depict the Cl- concentration (g/kgw) and velocity field after simulation time of 25, 50, and 100 years.
Table VI.3 Coordinates of sampling ports used to measure the Ca2+ concentration (reactive Henry problem simulation).
Table VI.4 Coordinates of sampling ports used to measure the seawater flux (reactive Henry problem simulation).
Figure.VI.7 Reactive Henry problem simulation using a kinetic coefficient kc=7.0×10-11 moles/s. 2D-porosity change (A), 2D-permeability change (B) and spatial distribution of the Cl- concentration (g/kgw) and velocity field (C) after 25, 50, and 100 years.
Fig. VI.7 displays the porosity and permeability change after 25, 50 and 100 years for kc= 7.0×10-11 moles/s. the spatial distribution of Cl- concentration and the veolocity field results of the development of porosity and permeability are shown in Fig. VI.7C. Maximum prosity increase rate is identified in the mixing zone and is about 1.07% after 100 years.
Figure.VI.8 Reactive Henry problem simulation using a kinetic coefficient kc=8.0×10-11 moles/s. 2D-porosity change (A), 2D-permeability change (B) and spatial distribution of the Cl- concentration (g/kgw) and velocity field (C) after 25, 50, and 100 years.
Fig. VI.8 displays the porosity and permeability change after 25, 50 and 100 years for kc= 8.0×10-11 moles/s. the spatial distribution of Cl- concentration and the veolocity field results of the development of porosity and permeability are shown in Fig. VI.8C. Maximum prosity increase rate is identified in the mixing zone and is about 1.22% after 100 years.
Figure.VI.9 Reactive Henry problem simulation using a kinetic coefficient kc=9.0×10-11 moles/s. 2D-porosity change (A), 2D-permeability change (B) and spatial distribution of the Cl- concentration (g/kgw) and velocity field (C) after 25, 50, and 100 years.
For kc=9×10-11 moles/s maximum porosity change is about 1.38% and maximum permeability change is about 2.5×10-10 m2 and located near to the discharge area (Fig. VI.9).
Figure.VI.10 Reactive Henry problem simulation using a kinetic coefficient kc=1.0×10-10 moles/s. 2D-porosity change (A), 2D-permeability change (B) and spatial distribution of the Cl- concentration (g/kgw) and velocity field (C) after 25, 50, and 100 years.
Figure.VI.11 Reactive Henry problem simulation using a kinetic coefficient kc=1.5×10-10 moles/s. 2D-porosity change (A), 2D-permeability change (B) and spatial distribution of the Cl- concentration (g/kgw) and velocity field (C) after 25, 50, and 100 years.
Figure.VI.12 Reactive Henry problem simulation using a kinetic coefficient kc=2.0×10-10 moles/s. 2D-porosity change (A), 2D-permeability change (B) and spatial distribution of the Cl- concentration (g/kgw) and velocity field (C) after 25, 50, and 100 years.
Fig. VI.10, Fig. VI.11 and Fig. VI.12 display the porosity and permeability development and their effect on the flow and transport for kc=10-10 moles/s, kc=1.5×10-9 moles/s and for kc=2.0×10-10 moles/s, respectively. Compared to the previous simulations, for a relatively high kc value the mixing zone is extended to the freshwater side and a displacement of the saltwater front to the freshwater side is identified. As the simulation progresses, there is a development of the sptatial distribution of Cl- concentration. Consequently, a movement of the saltwater front is created due to the development of the permeability in mixing zone which enhance more saltwater flux to the freshwater side and a larger mixing zone is created.
Figure.VI.13 Reactive Henry problem simulation using a kinetic coefficient kc=5.0×10-10 moles/s. Spatial distribution of the Cl- concentration (g/kgw) and velocity field after 25 (A), 50 (B), and 100 years (C).
Referring to the study of (Laabidi and Bouhlila 2015) and (Abarca et al. 2007b) we focus in this study on the following values:
Penetration of the seawater intrusion wedge noted Ltoe, as shown in Fig.VI.6 is calculated as the distance between the seaside boundary and the point where the 50% mixing isoline intersects the aquifer bottom.
qs is the seawater flux that enters the system through the seaside boundary (m3/m/s) integrated over the right boundary of the domain.
The result presented in table VI.5 shows the increase of Ltoe during time due to the development of velocity field during porosity and permeability change. Compared to the base case were Ltoe is about 0.6m the Ltoe in the reactive Henry problem can reach a values ranging from 0.965 to 1.881 m after 100 years for different kc values as presented in Table VI.5. On the other hand we can easily identify an increase of the seawater flux. Fig.VI.15 depicts an increase of the seawater flux qs to compare to the base case in sampling ports at selected points on the seaside boundary (Table VI.4). It can be seen from Fig. 15 that the model predicts reasonably sweater flux. Maximum value is reached in port N°4 and is about 1×51 10-5 m3/m/s for kc=2×10-10 moles/s where the hydraulic gradient is larger than in the other part of the seaside boundary, this can be caused by the increase of permeability and the development of the velocity field that enhance more seawater penetration to the freshwater side.
Table VI.5. The penetration of seawater (Ltoe) for different kc values and after simulation time of 25, 50 and 100 years.
To conclude, calcite dissolution induces an increase in the porosity which is related to the permeability. This would enhance further seawater flow and mixing in the simulated section. Chemical dissolution in such environment is another characteristic of real aquifers which the Henry problem does not account for. Calcite dissolution may affect the horizontal permeability which affects both seawater penetrations and the flux of saltwater that enters the aquifer through the seaside boundary. The shape and penetration of the saltwater intrusion wedge are controlled by the permeability which is related to the quantity of the dissolved calcite as presented above. However, it is difficult to evaluate the relative importance of this parameter. Here, we examine the results obtained when varying this parameter with respect to the base-case. For this reason several simulations were carried out to identify the role of the calcite dissolution using different kc values. To summarize, the penetration length increases and is mainly controlled by the horizontal permeability. The amount of dissolved calcite in the mixing zone increases with time inducing an increase in porosity due to the pore volume released by calcite dissolution phenomenon. The new hydrodynamic distribution controls the variation of the velocity field and the shape of the concentration repartition as shown in Fig. 13.
Figure.VI.14 Reactive Henry problem simulation: total aqueous calcium concentrations at locations of the sampling ports for different kinetic coefficient values.
Figure.VI.15 Reactive Henry problem simulation: Seawater flux at locations of the sampling ports for different kinetic coefficient values.
VI.II Alternative modeling approach
One of the most challenging in hydrogeochemical problems is the determination of analytical solution also called “closed solution’’ of reactive transport process in porous media. The complexity of the problem (irregular geometry, non-linearity of the mathematical formulation, coupling between two or more processes and complex boundary conditions) inhibited the development of analytical formalism.
Multicomponent reactive transport in porous media is modeled using numerical computer codes, and a few analytical solutions are available in literature (Bauer et al. 2001; Lunn et al. 1996; McLaren 1970; Montas 2003; Quezada et al. 2004; Sun et al. 1999). These analytical solutions are specific to radioactive decay and biochemical reactions. However, they are not extended to hydrogeochemical reactions where it is not trivial to formulate the interdependence (coupling) between the transport and the chemical process. A few of analytical solutions for multicomponent reactive transport in porous media have been developed and are as follow:
(Peter C. Lichtner et al. 1986) developed an exact solution to the diffusion reaction equation for single and two species based on the assumption of continuum representation of porous media and mineral reactions are considered to be in equilibrium with the fluid phase. The exact solution was compared to the result of a finite difference simulation of diffusion reaction of SiO2 (aq) in a porous column consisting of quartz in addition to an inert rock matrix.
(Phillips 1991) presented a new formalism based on the gradient reaction concept where he focused on the dissolution or precipitation of calcium carbonate in mixing zone. This formalism is based on the assumption that dissolution rate is fast; in such a way the solution can be considered in equilibrium with respect to calcite. This analytical solution is proportional to the square of the salinity gradient and to the second derivate of the calcium equilibrium concentration with respect to salinity.
(De Simoni et al. 2005b) presented a new procedure (general form of Phillips’s formalism) to solve multicomponent reactive transport problems in porous media where homogenous or classical heterogeneous reactions occur (Rubin 1983b). The methodology consists of a protocol composed of four steps: (1) definition of mobile conservative components, (2) solving the transport equation for conservative species, (3) evaluation of the concentration of aqueous species and (4) substituting the latter into the transport equation to evaluate the reaction rate. The reaction rate is proportional to mixing rate where u the vector of conservative components concentration and D is the dispersion tensor.
(Hayek et al. 2012) developed an analytical solution for one, two, and three dimensional diffusive transport coupled with dissolution or precipitation reactions and porosity change. Their solution consists in using the method of simplest equation method (Kudryashov 2005b, 2005a). Their solutions are exact and do not contain any approximation. They used this analytical solution for verifying numerical solution obtained from the reactive transport code OPENGEOSYS-GEMS (Kosakowski and Watanabe 2014b; Shao et al. 2009a; Shao et al. 2009b) and the COMSOL multiphysic software (http ://www.comsol.com/products/4.2/).
Recently (Anis Younes and Fahs 2013) developed a new semi-analytical solution for benchmarking reactive density driven flow models by extending the Henry saltwater intrusion problem to include dissolution and degradation reactions. This solution is based on the Fourier series formulation and incorporate about 6195 terms.
(Romanov and Dreybrodt 2006) used the analytical solution of (Phillips 1991) to set up an alternative modeling approach to determinate the porosity change in the saltwater-freshwater mixing zone of coastal carbonate aquifers. As a first step, they solve the advection transport equation for salinity s using SEAWAT program (Guo and Langevin 2002b) and the second step consists in calculating the dissolution capacity in the mixing zone by estimating the calcium equilibrium concentration that is obtained as a function of salinity using PHREEQC code (Parkhurst and Apello 1999b) for this purpose. They applied this approach for a stationary salinity distribution. Nevertheless porosity change caused a change of permeability which causes variations of salinity distribution. In this regard such processes should be treated step by step for a non-stationary distribution of salinity.
In this part of this chapter, we start from the work of (Romanov and Dreybrodt 2006) and a non-stationary procedure of Phillips’s formalism is implemented in GEODENS (Bouhlila 1999b; Bouhlila and Laabidi 2008a) to calculate the change of porosity after calculating the salinity gradient. First, we compare the simulation model results with the experimental work published by (Singurindy et al. 2004b) for both dissolution and precipitation reactions. Finally, this approach is used to model the effect of calcite dissolution in the seawater intrusion in coastal carbonate homogeneous aquifer.
VI.II.1. Mixing of carbonate waters
Alterations in mineral composition, precipitation or dissolution can occur during mixing of different waters. One common example is related to the mixing of saltwater and freshwater in coastal carbonate aquifer (Fig.VI.16). (L. N. Plummer 1975b) and (Wigley and Plummer 1976b) have provided a clear analysis of the geochemical process during mixing of freshwater with seawater in carbonate coastal aquifer. The most important result that they reached is that mixing of freshwater saturated with respect to calcite with seawater that is supersaturated with respect to calcite can leads to a solution undersaturated with respect to calcite. The resultant mixture is strongly dependent on temperature, on the nonlinear dependence of activity coefficients on ionic strength, on the nonlinear variation of equilibrium constant with temperature, on ionic strength effect, and most importantly on the partial pressure of CO2 (pCO2). Similarly, (Bottrell et al. 1991) showed that mixing of two carbonate groundwaters which are both saturated with respect to calcite can leads to undersaturated conditions if the original groundwaters have different pCO2.
Figure.VI.16. Mixing of different interstitial waters
Referring to (L. N. Plummer 1975b) the two most important factors generating the precipitation or the dissolution of such solids are the amount of the seawater in the mixture and the pCO2. If the fraction of the seawater in the mixture zone is low or the pCO2 is high, then these two conditions can generate the dissolution of the calcite in the groundwater. A real example of carbonate dissolution is found in coastal aquifers such as the Yucatan Peninsula in Mexico (Back et al. 1979b). If the fraction of the seawater in the mixture zone is high and pCO2 is low then, the precipitation of the calcite in the groundwater can be generated. These factors could explain the absence of dissolution in the coastal aquifer of Mallorca in Spain (Price and Herman 1991b). The undersaturation of the mixture can enhance the dissolution of calcite and thus increasing the concentration of Ca and HCO3, and may affect porosity and permeability. In areas which the coastal aquifer is mainly carbonated, calcite dissolution may lead to karstification and high permeability (Back et al. 1979b).
Fluid mixing is implemented either by titrating one solution into another solution or by specifying a mixing ratio between two solutions at one reaction step. In reactive transport models fluid mixing occurs due to the physical process of mechanical dispersion and molecular diffusion. The degree of mixing is therefore determined by the coefficient of hydrodynamic dispersion. A useful index of mixing between freshwater and seawater is the salinity s or the mixing ratio m=s/s0 where s0 is salinity of saltwater and s that of the mixture. During the mixture and without any chemical reaction, the concentration of all species in the mixture is evaluated by the following equation (in this case we focus in the concentration of calcium):
(VI.1)
where Cf and Cs are the concentration of freshwater and saltwater of calcium, respectively. As shown in Fig.VI.17 the straight line represent the concentration of Ca2+ in the mixture in equilibrium with respect to calcite. If a reaction between the mixture and the carbonate matrix occurs, the concentration of reacted Ca2+ is calculated as follows:
(VI.2)
Figure.VI.17. Mixing of seawater and freshwater saturated with respect to calcite.
As shown in Fig.VI.17, when seawater and freshwater, each initially saturated with respect to calcite, are mixed, the resulting calcium concentration lies along the straight line joining Cs and Cf. Unless Ceq(s) is precisely linear in s, the mixture may be unsaturated or supersaturated over the whole range of the salinity s if the curvature of Ceq(s) is everywhere concave downward or upward (Fig.VI.17 curve 1 and 2) or unsaturated over parts of the range and supersaturated over others (Fig.VI.17 curve 3).
As discussed before Cmix is a linear function of s, therefore
(VI.3)
can be calculated by using a geochemical code such as PHREEQC (Parkhurst and Apello 1999b) as a function of m as will be discussed in the next section.
VI.II.2.Alternative modeling approach: mathematical formulation and protocol of the numerical procedure
VI.II.2.1. Governing equation
GEODENS (Bouhlila 1999b; Bouhlila and Laabidi 2008a) is a 2D finite element code, that describes the transient convective and diffusive-dispersive transport of specie I, whose concentration is CI[MM-1], by the following equation
. (VI.4)
where ui[LT-1] is the effective velocity in the porous media and is related to the Darcy’s velocity qi[LT-1] by the following expression:
(VI.5)
where kij[L2] is the intrinsic permeability, [L3L-3] is the porosity. [ML-3] and [ML-1T-1] are the fluid density and viscosity, respectively. These last two are dependent on the solute concentrations. is the source term which represents the contribution of all geochemical reactions into the mass flux and Dij[L2T-1] is the diffusion-dispersion term and is expressed as follows:
(VI.6)
where is the Kronecker symbol, and are the longitudinal and the transversal dispersivity respectively, d0[L2T-1] is the molecular diffusion coefficient, ui[LT-1] is the velocity vector component in the i direction, uj[LT-1] is the velocity vector component in the j direction, and u[LT-1] is the norm of the velocity vector.
The coupling between density and concentration is assumed to be a linear function of the form:
(VI.7)
where is the fluid density at a base concentration, C0 ; and is a constant coefficient of density variability.
Eq. 7 is applied only for a solution that has salinity less than the seawater one. In this case the relation between fluid density and solute concentration can be expressed as a linear function (Herbert et al. 1988b) . If the fluid has salinity greater than typical seawater, then these equations may not be valid. In that case, a different empirical relation between salt concentration and fluid density should be used (Guo and Langevin 2002b).
If we consider that the dissolution or the precipitation rates are sufficiently fast, the solution can be assumed to be in equilibrium with respect to a given mineral. In this case the amount QI of dissolved or precipitated material per unit volume and unit time is given by:
(VI.8)
is the concentration of the Ith species in solution [MM-1]. It’s clear that if dissolution or precipitation reaction occurs Eq.8 becomes highly coupled, because become a function of all other ions concentrations.
For the remaining ions that do not contribute in the dissolution or the precipitation reactions of calcite such as Na, Cl, K, SO4… are described by the same transport equation but without source term
(VI.9)
Using Eq.9 GEODENS (Bouhlila 1999b; Bouhlila and Laabidi 2008a) can calculate chlorine concentration representative of all other conservative species j of seawater. If seawater mixes with freshwater saturated with respect to calcite, only the calcium concentration of the mixture (which is function of salinity as presented in Eq.1) remains in Eq.8.
If the salinity and the chemical composition of seawater and freshwater are defined, then the calcium concentration of the mixture equilibrated with respect to calcite Ceqmix can be calculated using PHREEQC (Parkhurst and Apello 1999b) as described before. Regarding Eq.2, then from the set of transport reaction equations only one remains (Romanov and Dreybrodt 2006):
(VI.10)
As presented by (Phillips 1991) and by differentiating with respect to the spatial coordinates and time
(VI.11)
The change of porosity is expressed as a function of as follows:
where M is the molecular weight of calcite and ρ the density of limestone.
VI.II.2.2. Implementation of porosity-permeability relationship
The approach used in this work to setup porosity-permeability relationship into the model is described in Fig.VI.18 and has been used by (Tianfu Xu et al. 2006a). The porosity is updated based on the mineral volume released from or accumulated in the porous media during dissolution or precipitation reactions respectively. The physicochemical processes are written at the macroscopic scale in GEODENS (Bouhlila 1999b; Bouhlila and Laabidi 2008a). The pore volume is continuously modified by the dissolution-precipitation reactions of salts. At every time step the quantity of the dissolved calcite is quantified, the change of porosity is calculated and the permeability is updated in every node of the mesh (Tianfu Xu et al. 2006a). The reaction rate which is the second derivate of the calcium equilibrium concentration in the equation is calculated by using PHREEQC (Parkhurst and Apello 1999b). For the next time step the same protocol is repeated by using the updated porosity and permeability distributions. The porosity at any node of the mesh and at every time step Δt is updated as follows (Tianfu Xu et al. 2006a):
(VI.12)
Where VCaCO3 [L3M-1] is the molar volume of calcite, and RCaCO3 [M L-3T-1] is the rate of calcite precipitation (-) or dissolution (+).
To describe the relation between porosity and permeability, capillary tubes or plane cracks models have been used and can give relatively simple relationships between these two parameters such as the Koseny-Carman formula (De Marsily 1981a; Tianfu Xu et al. 2006a) . (Sanford and Konikow 1989b) used an empirical relationship based on work done in reservoir engineering to represent the permeability variations due to the dissolution of calcite in a coastal aquifer as follows:
(VI.13)
According to these authors and for sedimentary rocks, a is about 20 and b is deduced from the initial known values of k and.
VI.II.2.3 Protocol used for the calculation
The protocol used in this work is based on the mixing function in PHREEQC (Parkhurst and Apello 1999b) as used by (Romanov and Dreybrodt 2006) extended to transient state and is defined as follows:
(a) Define the chemical composition of saline and freshwater solutions and equilibrate each of them with respect to a given partial pressure pCO2 of carbon dioxide.
(b) Mix the solutions with mixing ratio m of saline water.
(c) Equilibrate the mixture with respect to calcite and CO2 and find as a function of s.
(d) Find by use of Eq. 1 and Eq. 2.
(e) Fit the data points to a polynomial and calculate the second derivative by analytical differentiation with respect to s.
(f) Use the result in (e) to calculate the porosity and permeability change as presented in the flow chart in Fig.VI.18.
(g) Update the porosity and the permeability distribution to solve the flow and transport equation for the next time step.
Figure.VI.18. Flow chart of the numerical procedure for the alternative approach
VI.II.3. Testing the approach
(Singurindy et al. 2004b) conducted a laboratory experiment to calculate the total porosity change in a 2D-cell. The basic domain is 16 cm×16 cm×1 cm and two reservoirs supply it. Through the upside boundary freshwater in equilibrium with respect to calcite at pCO2 of 1 atm for the dissolution experiment and of 10-3.5 atm for the precipitation experiment is injected from the first reservoir under fixed head boundary of 1 cm. The second reservoir supplies the domain through the left hand side by a solution of water containing NaCl with concentration of 40 g/l for the dissolution experiment and of 80 g/l for the precipitation experiment under the same boundary conditions, in equilibrium with calcite and at pCO2 of 1 atm for the dissolution experiment and of 10-3.5 atm for the precipitation experiment. The two remaining boundaries present the outlets of the domain under fixed head boundary conditions equal to zero (h=0m). A dispersivity of 6.5 x 10-4 m and a porosity of 0.33 are used.
A grid of 441 elements is used to simulate this laboratory experiment. Fig.VI.19a and Fig.VI.19b show the calculated concentration distribution of the chloride (Cl-) in the modeled domain after the mixing of the two solutions for dissolution experiment and precipitation experiment respectively.
Figure.VI.19. Cl- isolines illustrate the calculated salinity distribution for (a) dissolution experiment and (b) precipitation experiment
Fig.VI.20a and Fig.VI.20b show the calcium concentration Ceq calculated by using PHREEQC (Parkhurst and Apello 1999b) for the dissolution and the precipitation experiment respectively. Fig.VI.21a and Fig.VI.21b shows the spatial distribution of the porosity change for the dissolution and the precipitation experiment respectively. We have found the same structure and location of the porosity change as in the simulation of (Romanov and Dreybrodt 2006) and the laboratory experiment of (Singurindy et al. 2004b). The laboratory experiment shows a porosity change along the domain diagonal (the mixing zone) as in the present work.
Figure.VI.20. Solubility of calcium carbonate as a function of salinity for (a) dissolution experiment and (b) precipitation experiment
As described before, the first simulation of the experimental work was carried out by (Romanov and Dreybrodt 2006). However, they did not take into account the effect of the change of porosity on the permeability. Actually, the changes of porosity cause a change in the permeability and this can be described by some relations, e.g. the Kozeny–Carman relation or the Sanford equation (section 3.2). In our work and in order to take into account the permeability evolution, for every time step we need to update the porosity and the permeability and recalculate the salinity distribution. This protocol is implemented in the current version of GEODENS (Bouhlila 1999b; Bouhlila and Laabidi 2008a).
Figure.VI.21. Spatial distribution of porosity change for (a) dissolution experiment and (b) precipitation experiment
The change of porosity has been determined as 0.240×10-4 year-1 for dissolution and as -0.050 year-1 for precipitation in the laboratory experiment of (Singurindy et al. 2004b). In our work, an average of 0.260×10-4 year-1 was determined for the dissolution experiment and an average of -0.053 year-1 for the precipitation experiment.
The results show that the approach presented in this work is in good agreement with the experiment and the previous numerical studies.
VI.II.4. Simulation and discussion : Henry Reactive problem:
In this study, the same geometry and parameters as in Henry problem is used (H. R. Henry 1964b). The Henry saltwater intrusion problem concerns a vertical slice through an isotropic, homogeneous, confined aquifer. For the flow boundary, a freshwater constant flux is applied from the inland boundary and a hydrostatic pressure is distributed on the sea boundary. As shown in Fig.VI.22 it is not easy to represent the seaside for the transport boundary. In this work we introduce a mixed boundary condition in the seaside to avoid difficulties in the finite element solution near the outflow region of the aquifer (Segol et al. 1975b). (Segol et al. 1975b) replaced the original homogeneous Dirichlet boundary by a mixed Dirichlet-Neumann boundary. The upper 20 cm was a zero concentration flux exit boundary, and the lower 80 cm remained a constant concentration Dirichlet boundary.
Figure.VI.22. Henry problem: domain and boundary conditions,(H. R. Henry 1964b).
(H. R. Henry 1964b) solved the governing equations in a nondimensional form and presented a semianalytical solution for the steady state distribution of salt concentration and stream function. Henry problem has been widely used as a test case of density dependent groundwater flow models benchmarked (Abarca et al. 2007b; Simpson and Clement 2003; C. I. Voss and Souza 1987). The boundary conditions of Henry’s problem are defined in Fig.VI.22. No flow occurs across the top and bottom boundaries. A freshwater flux is injected along the left vertical boundary of the domain. The prescribed pressure along the right vertical boundary is set at hydrostatic seawater pressure. The soil parameters used in Henry’s problem are summarized in Table VI.2. In the original Henry problem, a pseudo-permanent regime is reached after 100 minutes (Bouhlila 1999b; C. I. Voss 1984b). If the solid matrix is calcareous, the dissolution processes are perpetual and no permanent regime is reached.
Referring to the study of (Abarca et al. 2007b) seawater intrusion studies are usually concerned with three different variable characteristics of seawater intrusion that are calculated as follow:
Penetration of the seawater intrusion wedge noted Ltoe, as shown in Fig.VI.22 is calculated as the distance between the seaside boundary and the point where the 50% mixing isoline intersects the aquifer bottom.
Width of the mixing zone noted W, as shown in Fig.VI.22 is calculated as the average of WMZ (which is equal to the distance between isolines of 25% and 75% concentration) between 0.2xLtoe and 0.8xLtoe.
qs is the seawater flux that enters the system through the seaside boundary (m3/m/s) integrated over the right boundary of the domain.
Figure. VI.23. Reactive Henry problem simulation, 2D-porosity change computed using the alternative approach after 25 (a), 50 (b) and 100 years (c).
.
Figure. VI.24. Reactive Henry problem simulation, spatial distribution of the Cl- concentration (g/Kgw) and velocity field computed using the alternative approach after 25 (a), 50 (b) and 100 years (c).
Beside these three characteristics presented above, in this work, we focus on the non-stationary porosity development. During calcite dissolution there is a porosity changes causing an increase in permeability that would lead to increased seawater flux and mixing, thus enhancing further dissolution.
Table VI.6 The penetration of seawater (Ltoe) and averaged width of the mixing zone (W) for different simulation time.
The results presented by (Abarca et al. 2007b) and by (Nick et al. 2013b) showed that Ltoe, W and qs are sensitive to the transversal and longitudinal dispersivity. Their results showed that toe penetration of saltwater decreases with increasing transversal dispersivity. While width of the mixing zone and saltwater flux become larger.
Table VI.7 Coordinates of the sampling ports.
In this simulation we focus only on the effect of calcite dissolution on the behavior of saltwater intrusion such as the effect on the porosity and permeability and the consequence of this latter on the flow and transport. As expected, calcite dissolution can increase porosity and permeability due to the volume of calcite released the porous media. Fig.VI.23 depicts the porosity change after a simulation time of 25, 50 and 100 years. High values of porosity is about 8.8 10-5 year-1 and is found at the beginning of the outflow boundary (in the transition zone between seawater and mixed freshwater near the seaside boundary). This is reasonable, because the rate of calcite dissolution is maximal in the region where the gradient of salinity is maximal. The distribution of salinity and velocity field due to the change of porosity and permeability are illustrated in Fig.VI.24 after a simulation time of 25, 50 and 100 years. The result presented in table VI.6 shows the increase of Ltoe and W during time due to the development of velocity field during porosity and permeability change. Compared to the base case were Ltoe is about 0.6m and W is about 0.150; for the reactive Henry problem a value of 1.350m and 0.276 is reached after 100 years for Ltoe and W respectively. On the other hand we can easily identify an increase of the seawater flux. Fig.VI.25 depicts an increase of the seawater flux qs to compare to the base case in sampling ports at selected points on the seaside boundary (Table VI.7). It can be seen from Fig.VI.25 that the model predicts reasonably seawater flux, maximum value is reached in port N°4 and is about 1.72 10-5 m3/m/s where the hydraulic gradient is larger than in the other part of the seaside boundary, this can be caused by the increase of permeability and the development of the velocity field that enhance more seawater penetration to the freshwater side.
Figure.VI.25 Seawater flux at locations of the sampling ports.
To summarize, the penetration length increases and is mainly controlled by the horizontal permeability. The amount of calcite dissolved in the mixing zone increases with time inducing an increase in porosity due to the pore volume released by calcite dissolution reaction. The new hydrodynamic distribution controls the variation of the velocity field and the shape of the concentration repartition as shown in Fig.VI.24.
CHAPTER 7. CONCLUSIONS AND PERSPECTIVES
Results presented in this work show the interplay between geochemical reaction and transport in the mixing zone in coastal carbonate aquifers. The seawater intrusion phenomenon is a real risk on the subsurface water resources in several areas around the world. This risk can be increased by the solid matrix dissolution and precipitation processes in the saltwater–freshwater mixing zone in the case of limestone aquifer. Calcite dissolution and precipitation in the saltwater–freshwater mixing zone in coastal carbonate aquifers is described by a density dependent flow and multi species reactive transport model. This requires a set of highly nonlinearly coupled reaction–diffusion–advection equations. GEODENS code, used in this work can solve these equations by a finite element procedure. In order to quantify this risk, a relative complete modeling scheme is presented in this work. A finite element procedure is developed to solve these coupled equations and the geochemical processes are written at each element using Pitzer model for species activities calculation. At each time step, the geochemical fluxes are calculated with the previous time step results according to an explicit procedure. The porosity is calculated from the pore volume released or filled by calcite dissolution or precipitation in every node and at each time step. The new porosity value is used to update the permeability using an empirical equation. The new hydrodynamic distribution is the new inputs for the next time step to solve the problem in every nodes of the mesh.
In order to validate this model, different formulations and numerical schemes were used in the convection cell. The Galerkin method combined with the convective form of the transport equation CFTE using an implicit method provides results that compare very well with the results of Diersch (1988). The weighted residual approximation method developed here to evaluate the velocity is also satisfactory.
The salts and brine geochemistry model presented in this work is described according to the Pitzer model (Pitzer et al., 1984). It allows calculating the ion activities, the density of the solution, the factors concentrations of salts in the brine, and calculate the amounts of the dissolved or precipitated salts over time during dissolution-precipitation reactions which are described by a kinetic or equilibrium law. In order to validate the geochemical model, we have simulated the porosity change in a 1D diffusive water mixing problem (Rezeia et al.,2005) and we have showed that the results of the GEODENS code are satisfactory.
The hydrogeochemical model is the coupling between the density-dependent flow and transport model extended to ten species with the geochemical model where we retain a sequential iterations method between the physical and the chemical processes. In order to validate the coupled model, we have performed a simulation to predict the porosity evolution in a 2D flow cell. This result was compared to the experimental work presented by Singurindy and Berkowitz (2004).
Through the simulated example, we have shown that our hydrogeochemical model is in good agreement with the previous experimental work (Singurindy and Berkowitz, 2004) and the numerical one (Romanov and Dreybrodt, 2006).
VII.I Mixing inducing calcite dissolution
VII.I.1. Using the coupled hydrogeochemical model
Through the results of the simulations presented in this work, we have shown that the change of porosity is created in narrow regions at the freshwater side of the mixing zone. The behavior of the solution is discussed in terms of two parameters: velocity field which has a direct effect on the seawater flux, and the penetration length. An increase of the penetration length of the seawater is identified during calcite dissolution.
Sensitivity analyses were carried out in order to provide a clear analysis of the effect of the calcite dissolution on the flow and transport in porous media. As a first result, during calcite dissolution an increase of seawater flux is identified and is calculated in some sampling ports in the seaside boundary. Maximum seawater flux is located in region where seawater hydraulic head is larger than in the other part of the seaside boundary. Compared to the seawater flux calculated in the base case (nonreactive case), the increase of seawater flux from the seaside to the freshwater side is essentially the result of the development of the velocity field during permeability increase.
Dissolution processes in the mixing zone are depending on a large number of factors (mainly the differences in species concentrations, pCO2 and ionic strength of the boundary solution). Therefore, through the simulations carried out in this work the amount of calcite dissolution and consequently the porosity development are basically sensitive to the kinetic coefficient value. On the other hand, the largest dissolution rate is found mainly in two positions: one near the aquifer bottom and the other near the discharge area (the seaside boundary) in the mixing zone.
Another important result regarding the effect of calcite dissolution on the transport is the total aqueous calcium concentrations at locations of the sampling ports for different kinetic coefficient values. Generally during calcite dissolution an increase of Ca2+ concentration is expected. This can be seen clearly by comparing the measured concentration at sampling ports in the reactive simulation to the corresponding values detected for the base case (nonreactive case).
We essentially conclude that, for a real study case of coastal limestone aquifer the calcite dissolution processes are important to quantify the aquifer contaminated zone wedge and also its time evolution. On the other hand, it is highly recommended to make sure that the kinetic law of calcite dissolution is written carefully in a coupled model.
VII.I.2. Using anlutical solution
If dissolution or precipitation of calcite is sufficiently fast, such that the solution can be considered to be in equilibrium with respect to the dissolved or precipitated calcite and the equations describing the density dependent flow and reactive transport can be decoupled. In this regard, Phillips (1991) presented an analytical solution to solve such problems. This equation gives a better understanding of the location of dissolved or precipitated calcite during mixing of freshwater and seawater. Through the examples simulated in this work, we have shown that the approach gives reasonable results in comparison with the previous experiment work (Singurindy and Berkowitz, 2004) and numerical one (Romanov and Dreybrodt, 2006) both for dissolution and precipitation reactions. The result shows also that the changes of porosity are created in narrow regions at the freshwater side of the mixing zone.
VII.II Mixing inducing calcite precipitation
In this simulation, we have shown that during parallel flow deposition start around the inlet boundary and precipitation continue along the centerline cell (the mixing line). The effect of calcite deposition on the flow and transport is discussed in term of two parameters: porosity decrease and the total aqueous calcium concentrations in the same selected sampling ports. The code appears to capture the majority of the processes that can be identified during precipitation reaction such as the Ca2+ concentration, the amount of precipitated calcite and the porosity change but with different magnitude. These two parameters are the direct results of the geochemical phenomenon during the mixing process.
VII.III Perspectives
Further works are to focus on the following axis:
It would be interesting to integrate the thermal effects and their impact on reactive transport processes,
It would be of interest to perform laboratory experiments for calcite dissolution and precipitation reaction during mixing of seawater and freshwater and their impacts on the flow and transport
Integrating the 3rd dimension in GEODENS.
Modeling the aquifer heterogeneity. Aquifer heterogeneity strongly influences solute transport modeling.
APPENDIX
A.III.1. Integration of the fluid mass conservation equation
En utilisant l’approximation nodale de p et C, nous obtenons la forme matricielle suivante :
(II.15)
Avec :
A.III.2. Integration of the fluid mass conservation equation
The convective form of the transport equation CFTE
The diffusive form of the transport equation, DFTE:
A.III.3. Temporal Discretization
The discretization scheme is implicit by using a finite difference method (Bouhlila, 1999).
(II.19)
Si : Implicit Euler method méthode d’Euler implicite ;
: Semi-implicit Euler method;
: Explicit Euler method.
are the values of main variables fields (p and C) at t and respectively.
The two above matrix system founded in the integration of fluid mass conservation equation and transport equation can be written at as follow: (II.20)
An iterative procedure is adopted to solve the above nonlinear system.
A.III.4. Pitzer Coefficients values
Table A.III.1 Pitzer coefficients values (Michard 1989)
Table A.III.2 Pitzer coefficients values (Michard 1989)
Table A.III.3 Pitzer coefficients values (Michard 1989)
Table A.III.1 Pitzer coefficients values for density calculation (Monnin 1989)
REFERENCES
Abarca, E., Carrera, J., Sánchez-Vila, X., & Dentz, M. (2007b). Anisotropic dispersive Henry problem. Advances in Water Resources, 30(4), 913-926, doi:http://dx.doi.org/10.1016/j.advwatres.2006.08.005.
Abd-Elhamid, H. F., & Javadi, A. A. (2011). A density-dependant finite element model for analysis of saltwater intrusion in coastal aquifers. Journal of Hydrology, 401(3–4), 259-271, doi:http://dx.doi.org/10.1016/j.jhydrol.2011.02.028.
Ackerer, P., Younes, A., & Mose, R. (1999). Modeling Variable Density Flow and Solute Transport in Porous Medium: 1. Numerical Model and Verification. Transport in Porous Media, 35(3), 345-373, doi:10.1023/A:1006564309167.
Aguilera, D., Jourabchi, P., Spiteri, C., & Regnier, P. (2005). A knowledge‐based reactive transport approach for the simulation of biogeochemical dynamics in Earth systems. Geochemistry, Geophysics, Geosystems, 6(7).
Albets-Chico, X., & Kassinos, S. (2013). A consistent velocity approximation for variable-density flow and transport in porous media. Journal of Hydrology, 507(0), 33-51, doi:http://dx.doi.org/10.1016/j.jhydrol.2013.10.009.
Allison, J. D., Brown, D. S., & Novo-Gradac, K. J. (1991). MINTEQA2/ PRODEFA2, a geochemical assessment model for environmental systems, version 3.0 user's manual. . US Environmental Protectiona Agency Report EPA/600/3-91/021.
Appelo, C. A. J., Beekman, H. E., & Oosterbaan, A. W. A. (1984). Hydrochemistry of springs from dolomite reefs in the southern Alps of northern Italy. Int. Assoc. Hydrol.Sci. Publ., 150(1), 125-138, doi:10.2475/ajs.282.1.45.
Appelo, C. A. J., & Willemsen, A. (1987). Geochemical calculations and observations on salt water intrusions, I. A combined geochemical/minxing cell model. Journal of Hydrology, 94(3–4), 313-330, doi:http://dx.doi.org/10.1016/0022-1694(87)90058-8.
Arvidson, R. S., Ertan, I. E., Amonette, J. E., & Luttge, A. (2003). Variation in calcite dissolution rates:: A fundamental problem? Geochimica et Cosmochimica Acta, 67(9), 1623-1634, doi:http://dx.doi.org/10.1016/S0016-7037(02)01177-8.
Ataie-Ashtiani, B., Hassanizadeh, S. M., Oostrom, M., Celia, M. A., & White, M. D. (2001). Effective parameters for two-phase flow in a porous medium with periodic heterogeneities. Journal of Contaminant Hydrology, 49(1–2), 87-109, doi:http://dx.doi.org/10.1016/S0169-7722(00)00190-X.
Audigane, P., Gaus, I., Czernichowski-Lauriol, I., Pruess, K., & Xu, T. (2007). Two-dimensional reactive transport modeling of CO2 injection in a saline aquifer at the Sleipner site, North Sea. American Journal of Science, 307(7), 974-1008, doi:10.2475/07.2007.02.
Back, W., Hanshaw, B. B., Pyle, T. E., Plummer, L. N., & Weidie, A. E. (1979a). Geochemical significance of groundwater discharge and carbonate solution to the formation of Caleta Xel Ha, Quintana Roo, Mexico. Water Resources Research, 15(6), 1521-1535, doi:10.1029/WR015i006p01521.
Back, W., Hanshaw, B. B., Pyle, T. E., Plummer, L. N., & Weidie, A. E. (1979b). Geochemical significance of groundwater discharge and carbonate solution to the formation of Caleta Xel Ha, Quintana Roo, Mexico. Water Resources Research, 15(6), 1521-1535, doi:10.1029/WR015i006p01521.
Bakker, M., Oude Essink, G. H. P., & Langevin, C. D. (2004). The rotating movement of three immiscible fluids—a benchmark problem. Journal of Hydrology, 287(1–4), 270-278, doi:http://dx.doi.org/10.1016/j.jhydrol.2003.10.007.
Barlow, P., & Reichard, E. (2010). Saltwater intrusion in coastal regions of North America. Hydrogeology Journal, 18(1), 247-260, doi:10.1007/s10040-009-0514-3.
Bauer, P., Attinger, S., & Kinzelbach, W. (2001). Transport of a decay chain in homogenous porous media: analytical solutions. Journal of Contaminant Hydrology, 49(3–4), 217-239, doi:http://dx.doi.org/10.1016/S0169-7722(00)00195-9.
Bäverman, C., Strömberg, B., Moreno, L., & Neretnieks, I. (1999). CHEMFRONTS: a coupled geochemical and transport simulation tool. Journal of Contaminant Hydrology, 36(3–4), 333-351, doi:http://dx.doi.org/10.1016/S0169-7722(98)00152-1.
Bear, J. (1988). Dynamics of Fluids in Porous Media: Dover.
Benson, D. A., Carey, A. E., & Wheatcraft, S. W. (1998). Numerical advective flux in highly variable velocity fields exemplified by saltwater intrusion. Journal of Contaminant Hydrology, 34(3), 207-233, doi:http://dx.doi.org/10.1016/S0169-7722(98)00093-X.
Bethke, C. (1994). The Geochemist's Workbench. A User's Guide to Rxn, Act2, Tact, React and Gtplot. University of Illinois, Urbana-Champaign.
Bolton, E. W., Lasaga, A. C., & Rye, D. M. (1996). A model for the kinetic control of quartz dissolution and precipitation in porous media flow with spatially variable permeability: Formulation and examples of thermal convection. Journal of Geophysical Research: Solid Earth, 101(B10), 22157-22187, doi:10.1029/96JB01345.
Boluda-Botella, N., Gomis-Yagües, V., & Ruiz-Beviá, F. (2008). Influence of transport parameters and chemical properties of the sediment in experiments to measure reactive transport in seawater intrusion. Journal of Hydrology, 357(1–2), 29-41, doi:http://dx.doi.org/10.1016/j.jhydrol.2008.04.021.
Boris, S. K. (2001). Application of the Pitzer ion interaction model to natural hypersaline brines. Journal of Molecular Liquids, 91(1–3), 3-19, doi:http://dx.doi.org/10.1016/S0167-7322(01)00140-4.
Bottrell, S. H., Smart, P. L., Whitaker, F., & Raiswell, R. (1991). Geochemistry and isotope systematics of sulphur in the mixing zone of Bahamian blue holes. Applied Geochemistry, 6(1), 97-103, doi:http://dx.doi.org/10.1016/0883-2927(91)90066-X.
Boufadel, M. C. (2000). A mechanistic study of nonlinear solute transport in a groundwater–surface water system Under steady state and transient hydraulic conditions. Water Resources Research, 36(9), 2549-2565, doi:10.1029/2000WR900159.
Boufadel, M. C., Suidan, M. T., & Venosa, A. D. (1999). A numerical model for density-and-viscosity-dependent flows in two-dimensional variably saturated porous media. Journal of Contaminant Hydrology, 37(1–2), 1-20, doi:http://dx.doi.org/10.1016/S0169-7722(98)00164-8.
Boufadel, M. C., Xia, Y., & Li, H. (2011). Modeling solute transport and transient seepage in a laboratory beach under tidal influence. Environmental Modelling & Software, 26(7), 899-912, doi:http://dx.doi.org/10.1016/j.envsoft.2011.02.005.
Bouhlila, R. (1996). Ecoulements et transport couples dans les milieu poreux : modèles et applications. Les Annales Maghrébines de l’Ingénieur, 10(1–4), 7-29, doi:http://dx.doi.org/10.1016/j.jhydrol.2003.10.007.
Bouhlila, R. (1999a). Ecoulements, transports et réactions géochimiques couplés dans les milieux poreux. Cas des sels et des saumures. [Thèse de doctorat d'état es-sciences hydraulique]. Univ. Tunis El Manar.280 p., doi:http://10.13140/2.1.2270.0808 (in French).
Bouhlila, R. (1999b). Ecoulements, transports et réactions géochimiques couplés dans les milieux poreux. Cas des sels et des saumures. [Thèse de doctorat d'état es-sciences hydraulique]. 280, doi:http://10.13140/2.1.2270.0808 (in French).
Bouhlila, R., & Laabidi, E. (2008a). Impacts of calcite dissolution on seawater intrusion processes in coastal aquifers : density dependent flow and multi species reactive transport modelling. IAHS-AISH Publication, 320, 220-225.
Bouhlila, R., & Laabidi, E. (2008b). Impacts of calcite dissolution on seawater intrusion processes in coastal aquifers : density dependent flow and multi species reactive transport modelling. 320, V-356 p.
Buès, M. A., & Oltean, C. (2000). Numerical Simulations for Saltwater Intrusion by the Mixed Hybrid Finite Element Method and Discontinuous Finite Element Method. Transport in Porous Media, 40(2), 171-200, doi:10.1023/A:1006626230029.
Buhmann, D., & Dreybrodt, W. (1985). The kinetics of calcite dissolution and precipitation in geologically relevant situations of karst areas: 1. Open system. Chemical Geology, 48(1–4), 189-211, doi:http://dx.doi.org/10.1016/0009-2541(85)90046-4.
Busenberg, E., & Plummer, L. N. (1982). The kinetics of dissolution of dolomite in CO 2 -H 2 O systems at 1.5 to 65 degrees C and O to 1 atm PCO 2. American Journal of Science, 282(1), 45-78, doi:10.2475/ajs.282.1.45.
Caldeira, K., & Rau, G. H. (2000). Accelerating carbonate dissolution to sequester carbon dioxide in the ocean: Geochemical implications. Geophysical Research Letters, 27(2), 225-228, doi:10.1029/1999GL002364.
Caltagirone, J. P. (1975). Thermoconvective instabilities in a horizontal porous layer. Journal of Fluid Mechanics, 72(02), 269-287, doi:doi:10.1017/S0022112075003345.
Capaccioni, B., Didero, M., Paletta, C., & Salvadori, P. (2001). Hydrogeochemistry of groundwaters from carbonate formations with basal gypsiferous layers: an example from the Mt Catria–Mt Nerone ridge (Northern Appennines, Italy). Journal of Hydrology, 253(1–4), 14-26, doi:http://dx.doi.org/10.1016/S0022-1694(01)00480-2.
Carneiro, J., Boughriba, M., Correia, A., Zarhloule, Y., Rimi, A., & El Houadi, B. (2010). Evaluation of climate change effects in a coastal aquifer in Morocco using a density-dependent numerical model. Environmental Earth Sciences, 61(2), 241-252, doi:10.1007/s12665-009-0339-3.
Carrera, J., Vázquez-Suñé, E., Castillo, O., & Sánchez-Vila, X. (2004). A methodology to compute mixing ratios with uncertain end-members. Water Resources Research, 40(12), W12101, doi:10.1029/2003WR002263.
Centler, F., Shao, H., De Biase, C., Park, C.-H., Regnier, P., Kolditz, O., et al. (2010). GeoSysBRNS—A flexible multidimensional reactive transport model for simulating biogeochemical subsurface processes. Computers & Geosciences, 36(3), 397-405, doi:http://dx.doi.org/10.1016/j.cageo.2009.06.009.
Chen, J., Qian, H., & Li, P. (2013). Mixing Precipitation of CaCO3 in Natural Waters. Water, 5(4), 1712.
Cheng, A. H. D., & Ouazar, D. (2003). Coastal Aquifer Management-Monitoring, Modeling, and Case Studies: Taylor & Francis.
Cherubini, C., & Pastore, N. (2011). Critical stress scenarios for a coastal aquifer in southeastern Italy. Nat. Hazards Earth Syst. Sci., 11(5), 1381-1393, doi:10.5194/nhess-11-1381-2011.
Chu, M., Kitanidis, P. K., & McCarty, P. L. (2005). Modeling microbial reactions at the plume fringe subject to transverse mixing in porous media: When can the rates of microbial reaction be assumed to be instantaneous? Water Resources Research, 41(6), W06002, doi:10.1029/2004WR003495.
Cirpka, O. A., Frind, E. O., & Helmig, R. (1999). Numerical simulation of biodegradation controlled by transverse mixing. Journal of Contaminant Hydrology, 40(2), 159-182, doi:http://dx.doi.org/10.1016/S0169-7722(99)00044-3.
Cirpka, O. A., & Valocchi, A. J. (2007). Two-dimensional concentration distribution for mixing-controlled bioreactive transport in steady state. Advances in Water Resources, 30(6–7), 1668-1679, doi:http://dx.doi.org/10.1016/j.advwatres.2006.05.022.
Corbella, M., Ayora, C., & Cardellach, E. (2003). Dissolution of deep carbonate rocks by fluid mixing: a discussion based on reactive transport modeling. Journal of Geochemical Exploration, 78–79(0), 211-214, doi:http://dx.doi.org/10.1016/S0375-6742(03)00032-3.
Croucher, A. E., & O'Sullivan, M. J. (1995). The Henry Problem for Saltwater Intrusion. Water Resources Research, 31(7), 1809-1814, doi:10.1029/95WR00431.
Dagan, G., & Bear, J. (1968). Solving The Problem Of Local Interface Upconing In A Coastal Aquifer By The Method Of Small Perturbations. Journal of Hydraulic Research, 6(1), 15-44, doi:10.1080/00221686809500218.
Datta, B., Vennalakanti, H., & Dhar, A. (2009). Modeling and control of saltwater intrusion in a coastal aquifer of Andhra Pradesh, India. Journal of Hydro-environment Research, 3(3), 148-159, doi:http://dx.doi.org/10.1016/j.jher.2009.09.002.
De Marsily, G. (1981a). Quantitative hydrogeology: Ed. Masson (in French).
De Marsily, G. (1981b). Quantitative hydrogeology.: Ed. Masson (in French).
De Simoni, M., Carrera, J., Sánchez-Vila, X., & Guadagnini, A. (2005a). A procedure for the solution of multicomponent reactive transport problems. Water Resources Research, 41(11), W11410, doi:10.1029/2005WR004056.
De Simoni, M., Carrera, J., Sánchez-Vila, X., & Guadagnini, A. (2005b). A procedure for the solution of multicomponent reactive transport problems. Water Resources Research, 41(11), W11410, doi:10.1029/2005WR004056.
De Simoni, M., Sanchez-Vila, X., Carrera, J., & Saaltink, M. W. (2007). A mixing ratios-based formulation for multicomponent reactive transport. Water Resources Research, 43(7), W07419, doi:10.1029/2006WR005256.
Dhatt, G., & Touzot, G. (1981). Une présentation de la méthode des éléments finis: Maloine.
Diersch, H. J. (1988). Finite element modelling of recirculating density-driven saltwater intrusion processes in groundwater. Advances in Water Resources, 11(1), 25-43, doi:http://dx.doi.org/10.1016/0309-1708(88)90019-X.
Diersch, H. J., & Nillert., P. (1987). Saltwater intrusion processes in groundwater: novel computer simulations, field studies and interception techniques. Int. Symp. on Groundwater Monitoring and Management, 13(6), doi:10.1007/s10040-004-0333-5.
Diersch, H. J., Prochnow, D., & Thiele, M. (1984). Finite-element analysis of dispersion-affected saltwater upconing below a pumping well. Applied Mathematical Modelling, 8(5), 305-312, doi:http://dx.doi.org/10.1016/0307-904X(84)90143-4.
Diersch, H. J. G. (1996). Interactive, graphics-based finite-element simulation system FEFLOW for modeling groundwater flow, contaminant mass and heat transport processes, FEFLOW User's Manual Version 4.5. (Vol. 12, pp. 674-687): Institute for Water Resources Planning and System Research, Ltd.
Diersch, H. J. G., & Kolditz, O. (2002). Variable-density flow and transport in porous media: approaches and challenges. Advances in Water Resources, 25(8–12), 899-944, doi:http://dx.doi.org/10.1016/S0309-1708(02)00063-5.
Dreybrodt, W. (1980). Kinetics of the dissolution of calcite and its applications to karstification. Chemical Geology, 31(0), 245-269, doi:http://dx.doi.org/10.1016/0009-2541(80)90089-3.
El-Bihery, M. (2009). Groundwater flow modeling of Quaternary aquifer Ras Sudr, Egypt. Environmental Geology, 58(5), 1095-1105, doi:10.1007/s00254-008-1589-1.
Elder, J. W. (1967). Transient convection in a porous medium. Journal of Fluid Mechanics, 27(03), 609-623, doi:doi:10.1017/S0022112067000576.
Emmanuel, S., & Berkowitz, B. (2005). Mixing-induced precipitation and porosity evolution in porous media. Advances in Water Resources, 28(4), 337-344, doi:http://dx.doi.org/10.1016/j.advwatres.2004.11.010.
Engesgaard, P., & Kipp, K. L. (1992). A geochemical transport model for redox-controlled movement of mineral fronts in groundwater flow systems: A case of nitrate removal by oxidation of pyrite. Water Resources Research, 28(10), 2829-2843, doi:10.1029/92WR01264.
Engesgaard, P., & Traberg, R. (1996). Contaminant Transport at a Waste Residue Deposit: 2. Geochemical Transport Modeling. Water Resources Research, 32(4), 939-951, doi:10.1029/95WR03822.
Eriksson, N., & Destouni, G. (1997). Combined effects of dissolution kinetics, secondary mineral precipitation, and preferential flow on copper leaching from mining waste rock. Water Resources Research, 33(3), 471-483, doi:10.1029/96WR03466.
Fein, E. (1998). D3F – A simulator for density driven flow modelling, User's Manual. (Vol. 12, pp. 674-687): GRS. Braunscheig, Germany.
Felmy, A., Girvin, D., & Jenne, E. (1984). MINTEQ: A computer program for calculating aqueous geochemical equilibria.: U.S. Environmental Protection Agency, Washington, D.C., EPA/600/3-84/032.
Flukiger, F., & Bernard, D. (2009). A new numerical model for pore scale dissolution of calcite due to CO2 saturated water flow in 3D realistic geometry: Principles and first results. Chemical Geology, 265(1–2), 171-180, doi:http://dx.doi.org/10.1016/j.chemgeo.2009.05.004.
Freedman, V., & Ibaraki, M. (2002). Effects of chemical reactions on density-dependent fluid flow: on the numerical formulation and the development of instabilities. Advances in Water Resources, 25(4), 439-453, doi:http://dx.doi.org/10.1016/S0309-1708(01)00056-2.
Freedman, V. L., & Ibaraki, M. (2003). Coupled reactive mass transport and fluid flow: Issues in model verification. Advances in Water Resources, 26(1), 117-127, doi:http://dx.doi.org/10.1016/S0309-1708(02)00106-9.
Frind, E. O. (1982). Simulation of long-term transient density-dependent transport in groundwater. Advances in Water Resources, 5(2), 73-88, doi:http://dx.doi.org/10.1016/0309-1708(82)90049-5.
Galeati, G., & Gambolati, G. (1989). On boundary conditions and point sources in the finite element integration of the transport equation. Water Resources Research, 25(5), 847-856, doi:10.1029/WR025i005p00847.
Galeati, G., Gambolati, G., & Neuman, S. P. (1992). Coupled and partially coupled Eulerian-Lagrangian Model of freshwater-seawater mixing. Water Resources Research, 28(1), 149-165, doi:10.1029/91WR01927.
Gambolati, G., Putti, M., & Paniconi, C. (1999). Three-dimensional model of coupled density- dependent flow and miscible transport in groundwater. In: Bear et al (eds) Seawater intrusion in coastal aquifers:concepts, methods, and practices. (pp. 315-362): Kluwer, Dordrecht, The Netherlands.
Ghanbari, S., Al-Zaabi, Y., Pickup, G. E., Mackay, E., Gozalpour, F., & Todd, A. C. (2006). Simulation of CO2 Storage In Saline Aquifers. Chemical Engineering Research and Design, 84(9), 764-775, doi:http://dx.doi.org/10.1205/cherd06007.
Giambastiani, B. M. S., Antonellini, M., Oude Essink, G. H. P., & Stuurman, R. J. (2007). Saltwater intrusion in the unconfined coastal aquifer of Ravenna (Italy): A numerical model. Journal of Hydrology, 340(1–2), 91-104, doi:http://dx.doi.org/10.1016/j.jhydrol.2007.04.001.
Glassley, W. E., Nitao, J. J., & Grant, C. W. (2003). Three-dimensional spatial variability of chemical properties around a monitored waste emplacement tunnel. Journal of Contaminant Hydrology, 62–63(0), 495-507, doi:http://dx.doi.org/10.1016/S0169-7722(02)00153-5.
Gledhill, D. K., & Morse, J. W. (2006). Calcite solubility in Na–Ca–Mg–Cl brines. Chemical Geology, 233(3–4), 249-256, doi:http://dx.doi.org/10.1016/j.chemgeo.2006.03.006.
Gossel, W., Sefelnasr, A., & Wycisk, P. (2010). Modelling of paleo-saltwater intrusion in the northern part of the Nubian Aquifer System, Northeast Africa. Hydrogeology Journal, 18(6), 1447-1463, doi:10.1007/s10040-010-0597-x.
Goswami, R. R., & Clement, T. P. (2007). Laboratory-scale investigation of saltwater intrusion dynamics. Water Resources Research, 43(4), W04418, doi:10.1029/2006WR005151.
Graf, T., & Therrien, R. (2005). Variable-density groundwater flow and solute transport in porous media containing nonuniform discrete fractures. Advances in Water Resources, 28(12), 1351-1367, doi:http://dx.doi.org/10.1016/j.advwatres.2005.04.011.
Grasby, S. E., & Betcher, R. N. (2002). Regional hydrogeochemistry of the carbonate rock aquifer, southern Manitoba. Canadian Journal of Earth Sciences, 39(7), 1053-1063, doi:doi:10.1139/e02-021.
Greenberg, J. P., & Møller, N. (1989). The prediction of mineral solubilities in natural waters: A chemical equilibrium model for the Na-K-Ca-Cl-SO4-H2O system to high concentration from 0 to 250°C. Geochimica et Cosmochimica Acta, 53(10), 2503-2518, doi:http://dx.doi.org/10.1016/0016-7037(89)90124-5.
Grindrod, P., & Takase, H. (1996). Reactive chemical transport within engineered barriers. Journal of Contaminant Hydrology, 21(1–4), 283-296, doi:http://dx.doi.org/10.1016/0169-7722(95)00054-2.
Guo, W., & Langevin, C. D. (2002a). SEAWAT – User's Guide to SEAWAT: A computer program for simulation of three-dimensional variable-density ground-water flow. U.S.G.S. Techniques of Water-Resources Investigations, Book 6, Chapter A7, Tallahassee, Florida, 77pp, doi:10.1007/s10040-004-0333-5.
Guo, W., & Langevin, C. D. (2002b). SEAWAT – User's Guide to SEAWAT: A computer program for simulation of three-dimensional variable-density ground-water flow. U.S.G.S. Techniques of Water-Resources Investigations, Tallahassee, Florida, Open-File Report 01-434, 77pp.
Hanor, J. S. (1978). Precipitation of beachrock cements; mixing of marine and meteoric waters vs. CO 2 -degassing. Journal of Sedimentary Research, 48(2), 489-501, doi:10.1306/212f74b4-2b24-11d7-8648000102c1865d.
Hanshaw, B. B., & Back, W. (1979). Major geochemical processes in the evolution of carbonate—Aquifer systems. Journal of Hydrology, 43(1–4), 287-312, doi:http://dx.doi.org/10.1016/0022-1694(79)90177-X.
Harvie, C. E., Møller, N., & Weare, J. H. (1984). The prediction of mineral solubilities in natural waters: The Na-K-Mg-Ca-H-Cl-SO4-OH-HCO3-CO3-CO2-H2O system to high ionic strengths at 25°C. Geochimica et Cosmochimica Acta, 48(4), 723-751, doi:http://dx.doi.org/10.1016/0016-7037(84)90098-X.
Harvie, C. E., & Weare, J. H. (1980). The prediction of mineral solubilities in natural waters: the NaKMgCaClSO4H2O system from zero to high concentration at 25° C. Geochimica et Cosmochimica Acta, 44(7), 981-997, doi:http://dx.doi.org/10.1016/0016-7037(80)90287-2.
Hayek, M., Kosakowski, G., Jakob, A., & Churakov, S. V. (2012). A class of analytical solutions for multidimensional multispecies diffusive transport coupled with precipitation-dissolution reactions and porosity changes. Water Resources Research, 48(3), W03525, doi:10.1029/2011WR011663.
Hellevang, H., Aagaard, P., Oelkers, E. H., & Kvamme, B. (2005). Can Dawsonite Permanently Trap CO2? Environmental Science & Technology, 39(21), 8281-8287, doi:10.1021/es0504791.
Henderson, T. H., Mayer, K. U., Parker, B. L., & Al, T. A. (2009). Three-dimensional density-dependent flow and multicomponent reactive transport modeling of chlorinated solvent oxidation by potassium permanganate. Journal of Contaminant Hydrology, 106(3–4), 195-211, doi:http://dx.doi.org/10.1016/j.jconhyd.2009.02.009.
Henry, D., Touihri, R., Bouhlila, R., & Ben Hadid, H. (2012). Multiple flow solutions in buoyancy induced convection in a porous square box. Water Resources Research, 48(10), W10538, doi:10.1029/2012WR011995.
Henry, H. R. (1964a). Effect of dispersion on salt encroachment in coastal aquifers. US Geol Surv, Water-Supply Pap. No. 1613-C, 70–84, doi:10.1007/s10040-004-0333-5.
Henry, H. R. (1964b). Effect of dispersion on salt encroachment in coastal aquifers. US Geol Surv, Water-Supply Pap. No. 1613-C, 71–84.
Herbert, A. W., Jackson, C. P., & Lever, D. A. (1988a). Coupled groundwater flow and solute transport with fluid density strongly dependent upon concentration. Water Resources Research, 24(10), 1781-1795, doi:10.1029/WR024i010p01781.
Herbert, A. W., Jackson, C. P., & Lever, D. A. (1988b). Coupled groundwater flow and solute transport with fluid density strongly dependent upon concentration. Water Resources Research, 24(10), 1781-1795, doi:10.1029/WR024i010p01781.
Herman, J. S., & White, W. B. (1985). Dissolution kinetics of dolomite: Effects of lithology and fluid flow velocity. Geochimica et Cosmochimica Acta, 49(10), 2017-2026, doi:http://dx.doi.org/10.1016/0016-7037(85)90060-2.
Hess, J. W., & White, W. B. (1993). Groundwater geochemistry of the carbonate karst aquifer, southcentral Kentucky, U.S.A. Applied Geochemistry, 8(2), 189-204, doi:http://dx.doi.org/10.1016/0883-2927(93)90034-E.
Hickey, J. J. (1984). Field Testing the Hypothesis of Darcian Flow Through a Carbonate Aquifer. Ground Water, 22(5), 544-547, doi:10.1111/j.1745-6584.1984.tb01423.x.
Hilgers, C., & Urai, J. L. (2002). Experimental study of syntaxial vein growth during lateral fluid flow in transmitted light: first results. Journal of Structural Geology, 24(6–7), 1029-1043, doi:http://dx.doi.org/10.1016/S0191-8141(01)00089-X.
Holzbecher, E. (1998). Comments on “Constant-concentration boundary condition: Lessons from the HYDROCOIN variable-density groundwater benchmark problem” by L. F. Konikow, W. E. Sanford, and P. J. Campbell. Water Resources Research, 34(10), 2775-2778, doi:10.1029/98WR02190.
Holzbecher, E. (1998). Modeling density-driven flow in porous media. Principles, numerics, software, Springer Verlag, Berlin Heidelberg, 286pp, doi:10.1007/s10040-004-0333-5.
Hosokawa, T., Jinno, K., & Momii, K. (1990). Estimation of transverse dispersivity in the mixing zone of fresh-salt groundwater. Calibration and Reliability in Groundwater Modeling. IAHS Publ, 195, 149-158, doi:10.1007/s10040-004-0333-5.
Hui, Q., & Peiyue, L. (2011). Mixing Corrosion of CaCO3 in Natural Waters. E-Journal of Chemistry, 8(3), doi:10.1155/2011/891053.
Hui, Q., & Peiyue, L. (2012). Proportion dependent mixing effects of CaCO3 in natural waters. Asian Journal of Chemistry, 24(5), 2257.
Huisman, I. L. (1954). La formation des cones d'eau saumatre. IUGG IAHS-Publ., 37, 146-150, doi:10.1007/s10040-004-0333-5.
Huyakorn, P. S., Andersen, P. F., Mercer, J. W., & White, H. O. (1987). Saltwater intrusion in aquifers: Development and testing of a three-dimensional finite element model. Water Resources Research, 23(2), 293-312, doi:10.1029/WR023i002p00293.
Ibaraki, M. (1998). A robust and efficient numerical model for analyses of density-dependent flow in porous media. Journal of Contaminant Hydrology, 34(3), 235-246, doi:http://dx.doi.org/10.1016/S0169-7722(98)00092-8.
Jacobson, A. D., & Wasserburg, G. J. (2005). Anhydrite and the Sr isotope evolution of groundwater in a carbonate aquifer. Chemical Geology, 214(3–4), 331-350, doi:http://dx.doi.org/10.1016/j.chemgeo.2004.10.006.
Jakovovic, D., Werner, A. D., & Simmons, C. T. (2011). Numerical modelling of saltwater up-coning: Comparison with experimental laboratory observations. Journal of Hydrology, 402(3–4), 261-273, doi:http://dx.doi.org/10.1016/j.jhydrol.2011.03.021.
Johannsen, K., Kinzelbach, W., Oswald, S., & Wittum, G. (2002). The saltpool benchmark problem – numerical simulation of saltwater upconing in a porous medium. Advances in Water Resources, 25(3), 335-348, doi:http://dx.doi.org/10.1016/S0309-1708(01)00059-8.
Johns, R. T., & Rivera, A. (1996). Comment on “Dispersive Transport Dynamics in a Strongly Coupled Groundwater-Brine Flow System” by Curtis M. Oldenburg and Karsten Pruess. Water Resources Research, 32(11), 3405-3410, doi:10.1029/96WR02495.
Kang, Q., Lichtner, P., Viswanathan, H., & Abdel-Fattah, A. (2010). Pore Scale Modeling of Reactive Transport Involved in Geologic CO2 Sequestration. Transport in Porous Media, 82(1), 197-213, doi:10.1007/s11242-009-9443-9.
Katz, G. E., Berkowitz, B., Guadagnini, A., & Saaltink, M. W. (2011). Experimental and modeling investigation of multicomponent reactive transport in porous media. Journal of Contaminant Hydrology, 120–121(0), 27-44, doi:http://dx.doi.org/10.1016/j.jconhyd.2009.11.002.
Kerrou, J., & Renard, P. (2010). A numerical analysis of dimensionality and heterogeneity effects on advective dispersive seawater intrusion processes. Hydrogeology Journal, 18(1), 55-72, doi:10.1007/s10040-009-0533-0.
Kim, J., Schwartz, F. W., Xu, T., Choi, H., & Kim, I. S. (2004). Coupled Processes of Fluid Flow, Solute Transport, and Geochemical Reactions in Reactive Barriers. Vadose Zone J., 3(3), 867-874, doi:10.2136/vzj2004.0867.
Kipp, K. L., Jr. (1986). HST3D, A Computer Code for Simulation of Heat and Solute Transport in Three-dimensional Groundwater Flow Systems. US Geol Surv Water Resour Invest, Report 86-4095, doi:10.1023/A:1006504326005.
Knutson, C. E., Werth, C. J., & Valocchi, A. J. (2005). Pore-scale simulation of biomass growth along the transverse mixing zone of a model two-dimensional porous medium. Water Resources Research, 41(7), W07007, doi:10.1029/2004WR003459.
Kohfahl, C., Sprenger, C., Herrera, J. B., Meyer, H., Chacón, F. F., & Pekdeger, A. (2008). Recharge sources and hydrogeochemical evolution of groundwater in semiarid and karstic environments: A field study in the Granada Basin (Southern Spain). Applied Geochemistry, 23(4), 846-862, doi:http://dx.doi.org/10.1016/j.apgeochem.2007.09.009.
Kolditz, O., & Bauer, S. (2004). A process-oriented approach to computing multi-field problems in porous media. Journal of Hydroinformatics, 6, 225-244.
Kolditz, O., Bauer, S., Bilke, L., Böttcher, N., Delfs, J. O., Fischer, T., et al. (2012). OpenGeoSys: an open-source initiative for numerical simulation of thermo-hydro-mechanical/chemical (THM/C) processes in porous media. Environmental Earth Sciences, 67(2), 589-599, doi:10.1007/s12665-012-1546-x.
Kolditz, O., Ratke, R., Diersch, H.-J. G., & Zielke, W. (1998). Coupled groundwater flow and transport: 1. Verification of variable density flow and transport models. Advances in Water Resources, 21(1), 27-46, doi:http://dx.doi.org/10.1016/S0309-1708(96)00034-6.
Konikow, L. F., & Goode, D. J. (1996). A three-dimensional method-of-characteristics solute-transport model (MOC3D). U.S.G.S. Water-Resources Investigations Report 96-4267, 87pp, doi:10.1007/s10040-004-0333-5.
Konikow, L. F., Sanford, W. E., & Campbell, P. J. (1997). Constant-concentration boundary condition: Lessons from the HYDROCOIN variable-density groundwater benchmark problem. Water Resources Research, 33(10), 2253-2261, doi:10.1029/97WR01926.
Konz, M., Younes, A., Ackerer, P., Fahs, M., Huggenberger, P., & Zechner, E. (2009). Variable-density flow in heterogeneous porous media — Laboratory experiments and numerical simulations. Journal of Contaminant Hydrology, 108(3–4), 168-175, doi:http://dx.doi.org/10.1016/j.jconhyd.2009.07.005.
Kosakowski, G., & Watanabe, N. (2014a). OpenGeoSys-Gem: A numerical tool for calculating geochemical and porosity changes in saturated and partially saturated media. Physics and Chemistry of the Earth, Parts A/B/C, 70–71(0), 138-149, doi:http://dx.doi.org/10.1016/j.pce.2013.11.008.
Kosakowski, G., & Watanabe, N. (2014b). OpenGeoSys-Gem: A numerical tool for calculating geochemical and porosity changes in saturated and partially saturated media. Physics and Chemistry of the Earth, Parts A/B/C, 70–71(0), 138-149, doi:http://dx.doi.org/10.1016/j.pce.2013.11.008.
Kourakos, G., & Mantoglou, A. (2009). Pumping optimization of coastal aquifers based on evolutionary algorithms and surrogate modular neural network models. Advances in Water Resources, 32(4), 507-521, doi:http://dx.doi.org/10.1016/j.advwatres.2009.01.001.
Kourakos, G., & Mantoglou, A. (2011). Simulation and Multi-Objective Management of Coastal Aquifers in Semi-Arid Regions. Water Resources Management, 25(4), 1063-1074, doi:10.1007/s11269-010-9677-x.
Krause, R. E., & Clarke, J. S. (2004). Coastal ground water at risk – Saltwater contamination at Brunswick, Georgia and Hilton Head island, South Carolina. U.S. Geol. Surv.Water-Res. Invest. , Report 01-4107, doi:10.1007/s10040-004-0333-5.
Kudryashov, N. A. (2005a). Exact solitary waves of the Fisher equation. Physics Letters A, 342(1–2), 99-106, doi:http://dx.doi.org/10.1016/j.physleta.2005.05.025.
Kudryashov, N. A. (2005b). Simplest equation method to look for exact solutions of nonlinear differential equations. Chaos, Solitons & Fractals, 24(5), 1217-1231, doi:http://dx.doi.org/10.1016/j.chaos.2004.09.109.
Kulik, D., Wagner, T., Dmytrieva, S., Kosakowski, G., Hingerl, F., Chudnenko, K., et al. (2013). GEM-Selektor geochemical modeling package: revised algorithm and GEMS3K numerical kernel for coupled simulation codes. Computational Geosciences, 17(1), 1-24, doi:10.1007/s10596-012-9310-6.
Laabidi, E., & Bouhlila, R. (2015). Nonstationary porosity evolution in mixing zone in coastal carbonate aquifer using an alternative modeling approach. Environmental Science and Pollution Research, 22(13), 10070-10082, doi:10.1007/s11356-015-4207-2.
Laabidi, E., & Bouhlila, R. (2016a). Mixing inducing precipitation: impact of calcite precipitation reaction on the flow and transport. Carbonates and Evaporites (In Press).
Laabidi, E., & Bouhlila, R. (2016b). Reactive Henry problem: effect of calcite dissolution on seawater intrusion. Environmental Earth Sciences (In Press).
Laabidi, E., Gemma, B., & Elmer , V. d. B. (2010). Numerical model report Sarir-Tazerbo basins – Libya. Schlumberger Water Services Report (UK) Ltd The Pump House Coton Hill Shrewsbury SY1 2DP United Kingdom, doi:10.1007/s10040-004-0333-5.
Lagneau, V., Pipart, A., & Catalette, H. (2005). Modélisation couplée chimie-transport du comportement à long terme de la séquestration géologique de CO2 dans des aquifères salins profonds. Oil & Gas Science and Technology – Rev. IFP, 60(2), 231-247.
Langevin, C., Swain, E., & Wolfert, M. (2005). Simulation of integrated surface-water/ground-water flow and salinity for a coastal wetland and adjacent estuary. Journal of Hydrology, 314(1–4), 212-234, doi:http://dx.doi.org/10.1016/j.jhydrol.2005.04.015.
Langevin, C. D., Dausman, A. M., & Sukop, M. C. (2010). Solute and Heat Transport Model of the Henry and Hilleke Laboratory Experiment. Ground Water, 48(5), 757-770, doi:10.1111/j.1745-6584.2009.00596.x.
Le Gallo, Y., Bildstein, O., & Brosse, E. (1998). Coupled reaction-flow modeling of diagenetic changes in reservoir permeability, porosity and mineral compositions. Journal of Hydrology, 209(1–4), 366-388, doi:http://dx.doi.org/10.1016/S0022-1694(98)00183-8.
Lee, Y.-J., & Morse, J. W. (1999). Calcite precipitation in synthetic veins: implications for the time and fluid volume necessary for vein filling. Chemical Geology, 156(1–4), 151-170, doi:http://dx.doi.org/10.1016/S0009-2541(98)00183-1.
Legrand, H. E., & Stringfield, V. T. (1971). Development and Distribution of Permeability in Carbonate Aquifers. Water Resources Research, 7(5), 1284-1294, doi:10.1029/WR007i005p01284.
Li, H., & Boufadel, M. C. (2011). A tracer study in an Alaskan gravel beach and its implications on the persistence of the Exxon Valdez oil. Marine Pollution Bulletin, 62(6), 1261-1269, doi:http://dx.doi.org/10.1016/j.marpolbul.2011.03.011.
Li, H., Boufadel, M. C., & Weaver, J. W. (2008). Tide-induced seawater–groundwater circulation in shallow beach aquifers. Journal of Hydrology, 352(1–2), 211-224, doi:http://dx.doi.org/10.1016/j.jhydrol.2008.01.013.
Li, X., Hu, B. X., Burnett, W. C., Santos, I. R., & Chanton, J. P. (2009). Submarine Ground Water Discharge Driven by Tidal Pumping in a Heterogeneous Aquifer. Ground Water, 47(4), 558-568, doi:10.1111/j.1745-6584.2009.00563.x.
Lichtner, P. C. (2000). FLOTRAN: User's Manual, Los Alamos National Laboratory Report. (pp. 1689-1701): McPherson, B.J.O.L. and Chapman, D.S., 1996, Thermal Analysis of the Southern Powder River Basin, Wyoming, Geophysics.
Lichtner, P. C., & Felmy, A. R. (2003). Estimation of Hanford SX tank waste compositions from historically derived inventories. Computers & Geosciences, 29(3), 371-383, doi:http://dx.doi.org/10.1016/S0098-3004(03)00012-8.
Lichtner, P. C., Oelkers, E. H., & Helgeson, H. C. (1986). Exact and numerical solutions to the moving boundary problem resulting from reversible heterogeneous reactions and aqueous diffusion in a porous medium. Journal of Geophysical Research: Solid Earth, 91(B7), 7531-7544, doi:10.1029/JB091iB07p07531.
Lichtner, P. C., Yabusaki, S., Pruess, K., & Steefel, C. I. (2004). Role of Competitive Cation Exchange on Chromatographic Displacement of Cesium in the Vadose Zone beneath the Hanford S/SX Tank Farm. Vadose Zone J., 3(1), 203-219, doi:10.2136/vzj2004.2030.
Lin, J., Snodsmith, J. B., Zheng, C., & Wu, J. (2009). A modeling study of seawater intrusion in Alabama Gulf Coast, USA. Environmental Geology, 57(1), 119-130, doi:10.1007/s00254-008-1288-y.
Liu, C. W., & Narasimhan, T. N. (1989). Redox-controlled multiple-species reactive chemical transport: 2. Verification and application. Water Resources Research, 25(5), 883-910, doi:10.1029/WR025i005p00883.
Loucks, R. G. (1999). Paleocave carbonate reservoirs; origins, burial-depth modifications, spatial complexity, and reservoir implications. AAPG Bulletin, 83(11), 1795-1834.
Lu, C., Kitanidis, P. K., & Luo, J. (2009). Effects of kinetic mass transfer and transient flow conditions on widening mixing zones in coastal aquifers. Water Resources Research, 45(12), W12402, doi:10.1029/2008WR007643.
Lucia, F. J. (1999). Carbonate Reservoir Characterization: Springer.
Lunn, M., Lunn, R. J., & Mackayb, R. (1996). Determining analytic solutions of multiple species contaminant transport, with sorption and decay. Journal of Hydrology, 180(1–4), 195-210, doi:http://dx.doi.org/10.1016/0022-1694(95)02891-9.
Luyun Jr, R., Momii, K., & Nakagawa, K. (2009). Laboratory-scale saltwater behavior due to subsurface cutoff wall. Journal of Hydrology, 377(3–4), 227-236, doi:http://dx.doi.org/10.1016/j.jhydrol.2009.08.019.
MacQuarrie, K. T. B., & Mayer, K. U. (2005). Reactive transport modeling in fractured rock: A state-of-the-science review. Earth-Science Reviews, 72(3–4), 189-227, doi:http://dx.doi.org/10.1016/j.earscirev.2005.07.003.
Mangold, D. C., & Tsang, C.-F. (1991). A summary of subsurface hydrological and hydrochemical models. Reviews of Geophysics, 29(1), 51-79, doi:10.1029/90RG01715.
Mao, X., Enot, P., Barry, D. A., Li, L., Binley, A., & Jeng, D. S. (2006a). Tidal influence on behaviour of a coastal aquifer adjacent to a low-relief estuary. Journal of Hydrology, 327(1–2), 110-127, doi:http://dx.doi.org/10.1016/j.jhydrol.2005.11.030.
Mao, X., Prommer, H., Barry, D. A., Langevin, C. D., Panteleit, B., & Li, L. (2006b). Three-dimensional model for multi-component reactive transport with variable density groundwater flow. Environmental Modelling & Software, 21(5), 615-628, doi:http://dx.doi.org/10.1016/j.envsoft.2004.11.008.
Mayer, K. U., Frind, E. O., & Blowes, D. W. (2002). Multicomponent reactive transport modeling in variably saturated porous media using a generalized formulation for kinetically controlled reactions. Water Resources Research, 38(9), 1174, doi:10.1029/2001WR000862.
McIntosh, J. C., & Walter, L. M. (2006). Paleowaters in Silurian-Devonian carbonate aquifers: Geochemical evolution of groundwater in the Great Lakes region since the Late Pleistocene. Geochimica et Cosmochimica Acta, 70(10), 2454-2479, doi:http://dx.doi.org/10.1016/j.gca.2006.02.002.
McLaren, A. D. (1970). Temporal and vectorial reactions of nitrogen in soil: a review. Canadian Journal of Soil Science, 50(2), 97-109, doi:doi:10.4141/cjss70-017.
Michael, H. A., Mulligan, A. E., & Harvey, C. F. (2005). Seasonal oscillations in water exchange between aquifers and the coastal ocean. [10.1038/nature03935]. Nature, 436(7054), 1145-1148, doi:http://www.nature.com/nature/journal/v436/n7054/suppinfo/nature03935_S1.html.
Michard, G. (1989). Equilibres chimiques dans les eaux naturelles (Collection sciences et techniques ed.). Paris: Publisud.
Miller, C. W., & Benson, L. V. (1983). Simulation of solute transport in a chemically reactive heterogeneous system: Model development and application. Water Resources Research, 19(2), 381-391, doi:10.1029/WR019i002p00381.
Misut, P. E., & Voss, C. I. (2004). Simulation of seawater intrusion resulting from proposed expanded pumpage in New York City, USA. In T. M. Cass, & F. P. George (Eds.), Developments in Water Science (Vol. Volume 55, Part 2, pp. 1595-1606): Elsevier.
Monnin, C. (1989). An ion interaction model for the volumetric properties of natural waters: Density of the solution and partial molal volumes of electrolytes to high concentrations at 25°C. Geochimica et Cosmochimica Acta, 53(6), 1177-1188, doi:http://dx.doi.org/10.1016/0016-7037(89)90055-0.
Montas, H. J. (2003). An analytical solution of the three-component transport equation with application to third-order transport. Water Resources Research, 39(2), 1036, doi:10.1029/2002WR001288.
Moore, C. H. (1989). Carbonate Diagenesis and Porosity: Elsevier Science.
Moore, W. S. (2010). The Effect of Submarine Groundwater Discharge on the Ocean. Annual Review of Marine Science, 2(1), 59-88, doi:doi:10.1146/annurev-marine-120308-081019.
Morse, J. W., & Arvidson, R. S. (2002). The dissolution kinetics of major sedimentary carbonate minerals. Earth-Science Reviews, 58(1–2), 51-84, doi:http://dx.doi.org/10.1016/S0012-8252(01)00083-6.
Morse, J. W., & Berner, R. A. (1972). Dissolution kinetics of calcium carbonate in sea water; I, A kinetic origin for the lysocline. American Journal of Science, 272(9), 840-851, doi:10.2475/ajs.272.9.840.
Mulligan, A. E., Evans, R. L., & Lizarralde, D. (2007). The role of paleochannels in groundwater/seawater exchange. Journal of Hydrology, 335(3–4), 313-329, doi:http://dx.doi.org/10.1016/j.jhydrol.2006.11.025.
Nasri, N., Bouhlila, R., Saaltink, M. W., & Gamazo, P. (2015). Modeling the hydrogeochemical evolution of brine in saline systems: Case study of the Sabkha of Oum El Khialate in South East Tunisia. Applied Geochemistry, 55, 160-169, doi:http://dx.doi.org/10.1016/j.apgeochem.2014.11.003.
Nick, H. M., Raoof, A., Centler, F., Thullner, M., & Regnier, P. (2013a). Reactive dispersive contaminant transport in coastal aquifers: Numerical simulation of a reactive Henry problem. Journal of Contaminant Hydrology, 145(0), 90-104, doi:http://dx.doi.org/10.1016/j.jconhyd.2012.12.005.
Nick, H. M., Raoof, A., Centler, F., Thullner, M., & Regnier, P. (2013b). Reactive dispersive contaminant transport in coastal aquifers: Numerical simulation of a reactive Henry problem. Journal of Contaminant Hydrology, 145(0), 90-104, doi:http://dx.doi.org/10.1016/j.jconhyd.2012.12.005.
Nishikawa, T., Siade, A. J., Reichard, E. G., Ponti, D. J., Canales, A. G., & Johnson, T. A. (2009). Stratigraphic controls on seawater intrusion and implications for groundwater management, Dominguez Gap area of Los Angeles, California, USA. Hydrogeology Journal, 17(7), 1699-1725, doi:10.1007/s10040-009-0481-8.
Nurmi, R., & Standen, E. (1997). Carbonates, the inside story. Middle East Well Evaluation Review, 18, 27-41.
Oldenburg, C. M., & Pruess, K. (1995). Dispersive Transport Dynamics in a Strongly Coupled Groundwater-Brine Flow System. Water Resources Research, 31(2), 289-302, doi:10.1029/94WR02272.
Olivella, S., Gens, A., Carrera, J., & Alonso, E. E. (1996). Numerical formulation for a simulator (CODE_BRIGHT) for the coupled analysis of saline media. Engineering Computations, 13(7), 87-112, doi:doi:10.1108/02644409610151575.
Oswald, S. E., & Kinzelbach, W. (2004). Three-dimensional physical benchmark experiments to test variable-density flow models. Journal of Hydrology, 290(1–2), 22-42, doi:http://dx.doi.org/10.1016/j.jhydrol.2003.11.037.
Oude Essink, G. H. P. (2001). Improving fresh groundwater supply—problems and solutions. Ocean & Coastal Management, 44(5–6), 429-449, doi:http://dx.doi.org/10.1016/S0964-5691(01)00057-6.
Oude Essink, G. H. P. (2001). Saltwater intrusion in 3D large-scale aquifers: a dutch case. Physics and Chemistry of the Earth, Part B: Hydrology, Oceans and Atmosphere, 26(4), 337-344, doi:http://dx.doi.org/10.1016/S1464-1909(01)00016-8.
Oude Essink, G. H. P., van Baaren, E. S., & de Louw, P. G. B. (2010). Effects of climate change on coastal groundwater systems: A modeling study in the Netherlands. Water Resources Research, 46(10), W00F04, doi:10.1029/2009WR008719.
Oude Essink, G. P. (2001). Salt Water Intrusion in a Three-dimensional Groundwater System in The Netherlands: A Numerical Study. Transport in Porous Media, 43(1), 137-158, doi:10.1023/A:1010625913251.
Parkhurst, D. L., & Apello, C. A. J. (1999a). User’s Guide to PHREEQC – A Computer Program for Speciation, Reaction-path, 1D-transport, and Inverse Geochemical Calculations. In U. G. S. Technical Report 99-4259, USA (Ed.), (Version 2).
Parkhurst, D. L., & Apello, C. A. J. (1999b). User’s Guide to PHREEQC – A Computer Program for Speciation, Reaction-path, 1D-transport, and Inverse Geochemical Calculations. In U. G. S. Technical Report 99-4259, USA (Ed.), (Version 2): Technical Report 99-4259, US Geological Survey, USA.
Parkhurst, D. L., Kipp, K. L., Engesgaard, P., & Charlton, S. R. (2004). PHAST- a program for simulating ground-water flow and multicomponent geochemical reactions. (pp. 154): USGS Techniques and Methods 6-A8.
Pauline, H., Pascal, A., Julie, L., Philippe, N., & Vincent, L. (2011). Tracking and CO2 leakage from deep saline to fresh groundwaters: Development of sensitive monitoring techniques. Energy Procedia, 4(0), 3443-3449, doi:http://dx.doi.org/10.1016/j.egypro.2011.02.269.
Phillips, O. M. (1991). Geological fluid dynamics : sub-surface flow and reactions. Cambridge, UK ; New York: Cambridge University Press.
Pitzer, K. S. (1973). Thermodynamics of electrolytes. I. Theoretical basis and general equations. The Journal of Physical Chemistry, 77(2), 268-277, doi:10.1021/j100621a026.
Pitzer, K. S. (1979). Theory, ion interaction approach. In: Activity Coefficients in Electrolyte Solutions. Vol. 1. R.M. Pytkowicz (ed.). CRC Press,, 157-208.
Pitzer, K. S. (1991). On interaction approach: theory and data correlation. In 'Activity Coefficients in Electrolyte Solutions'. 2nd Ed., ed. K. S. Pitzer, CRC Press, Boca Raton, 279-434.
Pitzer, K. S., Christopher Peiper, J., & Busey, R. H. (1984). Thermodynamic Properties of Aqueous Sodium Chloride Solutions. Journal of Physical and Chemical Reference Data 13(1).
Plummer, L. N. (1975a). Mixing of Sea Water with Calcium Carbonate Ground Water. Geological Society of America Memoirs, 142, 219-236, doi:10.1130/MEM142-p219.
Plummer, L. N. (1975b). Mixing of sea water with calcium carbonate ground water. Geological Society of America Memoirs, 142, 219-236, doi:10.1130/MEM142-p219.
Plummer, L. N., Busby, J. F., Lee, R. W., & Hanshaw, B. B. (1990). Geochemical Modeling of the Madison Aquifer in Parts of Montana, Wyoming, and South Dakota. Water Resources Research, 26(9), 1981-2014, doi:10.1029/WR026i009p01981.
Plummer, L. N., Parkhurst, D. L., & Wigle, T. M. L. (1979). Critical Review of the Kinetics of Calcite Dissolution and Precipitation. In Chemical Modeling in Aqueous Systems (ACS Symposium Series ed., Vol. 93, pp. 537-573, ACS Symposium Series, Vol. 1–4).
Pool, M., & Carrera, J. (2010). Dynamics of negative hydraulic barriers to prevent seawater intrusion. Hydrogeology Journal, 18(1), 95-105, doi:10.1007/s10040-009-0516-1.
Post, V. E. A. (2005). Fresh and saline groundwater interaction in coastal aquifers: Is our technology ready for the problems ahead? Hydrogeology Journal, 13(1), 120-123, doi:10.1007/s10040-004-0417-2.
Poulet, T., Gross, L., Georgiev, D., & Cleverley, J. (2012). escript-RT: Reactive transport simulation in Python using escript. Computers & Geosciences, 45(0), 168-176, doi:http://dx.doi.org/10.1016/j.cageo.2011.11.005.
Praveena, S. M., & Aris, A. Z. (2010). Groundwater resources assessment using numerical model: A case study in low-lying coastal area. International Journal of Environmental Science & Technology, 7(1), 135-146, doi:10.1007/BF03326125.
Price, R. M., & Herman, J. S. (1991a). Geochemical investigation of salt-water intrusion into a coastal carbonate aquifer: Mallorca, Spain. Geological Society of America Bulletin, 103(10), 1270-1279, doi:10.1130/0016-7606(1991)103<1270:gioswi>2.3.co;2.
Price, R. M., & Herman, J. S. (1991b). Geochemical investigation of salt-water intrusion into a coastal carbonate aquifer: Mallorca, Spain. Geological Society of America Bulletin, 103(10), 1270-1279, doi:10.1130/0016-7606(1991)103<1270:gioswi>2.3.co;2.
Prommer, H. (2002). PHT3D – A Reactive Multicomponent Transport Model for Saturated Porous Media.Technical.: Contaminated Land Assessment and Remediation Research Centre, The University of Edinburgh.
Purser, B. H. (1980). Sédimentation et diagenèse des carbonates néritiques récents: les éléments de la sédimentation et de la diagenèse (Vol. vol. 1): Technip.
Putti, M., & Paniconi, C. (1995). Picard and Newton linearization for the coupled model for saltwater intrusion in aquifers. Advances in Water Resources, 18(3), 159-170, doi:http://dx.doi.org/10.1016/0309-1708(95)00006-5.
Quezada, C. R., Clement, T. P., & Lee, K.-K. (2004). Generalized solution to multi-dimensional multi-species transport equations coupled with a first-order reaction network involving distinct retardation factors. Advances in Water Resources, 27(5), 507-520, doi:http://dx.doi.org/10.1016/j.advwatres.2004.02.013.
Rechard, R. P. (1999). Historical relationship between performance assessment for radioactive waste disposal and other types of risk assessment. Risk Analysis, 19(5), 763-807, doi:10.1111/j.1539-6924.1999.tb00446.x.
Rege, S. D., & Fogler, H. S. (1989). Competition among flow, dissolution, and precipitation in porous media. AIChE Journal, 35(7), 1177-1185, doi:10.1002/aic.690350713.
Reilly, T. E., & Goodman, A. S. (1985). Quantitative analysis of saltwater-freshwater relationships in groundwater systems—A historical perspective. Journal of Hydrology, 80(1–2), 125-160, doi:http://dx.doi.org/10.1016/0022-1694(85)90078-2.
Reilly, T. E., & Goodman, A. S. (1987). Analysis of saltwater upconing beneath a pumping well. Journal of Hydrology, 89(3–4), 169-204, doi:http://dx.doi.org/10.1016/0022-1694(87)90179-X.
Rezaei, M., Sanz, E., Raeisi, E., Ayora, C., Vázquez-Suñé, E., & Carrera, J. (2005). Reactive transport modeling of calcite dissolution in the fresh-salt water mixing zone. Journal of Hydrology, 311(1–4), 282-298, doi:http://dx.doi.org/10.1016/j.jhydrol.2004.12.017.
Rifai, M. N. E., Kaufman, W. J., & Todd, D. K. (1956). Dispersion Phenomena in Laminar Flow Through Porous Media. Progress Report 2, Canal Seepage Rresearch. Univ. of California, Berkeley, doi:10.1007/s10040-004-0333-5.
Risacher, F., & Clement, A. (2001). A computer program for the simulation of evaporation of natural waters to high concentration. Computers & Geosciences, 27(2), 191-201, doi:http://dx.doi.org/10.1016/S0098-3004(00)00100-X.
Robinson, C., Li, L., & Barry, D. A. (2007). Effect of tidal forcing on a subterranean estuary. Advances in Water Resources, 30(4), 851-865, doi:http://dx.doi.org/10.1016/j.advwatres.2006.07.006.
Romanov, D., & Dreybrodt, W. (2006). Evolution of porosity in the saltwater–freshwater mixing zone of coastal carbonate aquifers: An alternative modelling approach. Journal of Hydrology, 329(3–4), 661-673, doi:http://dx.doi.org/10.1016/j.jhydrol.2006.03.030.
Rubin, J. (1983a). Transport of reacting solutes in porous media: Relation between mathematical nature of problem formulation and chemical nature of reactions. Water Resources Research, 19(5), 1231-1252, doi:10.1029/WR019i005p01231.
Rubin, J. (1983b). Transport of reacting solutes in porous media: Relation between mathematical nature of problem formulation and chemical nature of reactions. Water Resources Research, 19(5), 1231-1252, doi:10.1029/WR019i005p01231.
Rueedi, J., Purtschert, R., Beyerle, U., Alberich, C., & Kipfer, R. (2005). Estimating groundwater mixing ratios and their uncertainties using a statistical multi parameter approach. Journal of Hydrology, 305(1–4), 1-14, doi:http://dx.doi.org/10.1016/j.jhydrol.2004.06.044.
Saaltink, M. W., Ayora, C., Batlle, F., Carrera, J., & Olivella, S. (2004). RETRASO, a code for modeling reactive transport in saturated and unsaturated porous media. Geologica Acta, 2, 235-251.
Saaltink, M. W., Ayora, C., & Carrera, J. (1998). A mathematical formulation for reactive transport that eliminates mineral concentrations. Water Resources Research, 34(7), 1649-1656, doi:10.1029/98WR00552.
Saaltink, M. W., Carrera, J., & Ayora, C. (2001). On the behavior of approaches to simulate reactive transport. Journal of Contaminant Hydrology, 48(3–4), 213-235, doi:http://dx.doi.org/10.1016/S0169-7722(00)00172-8.
Saaltink, M. W., Pifarré, F. B., Ayora, C., Carrera, J., & Pastallé, S. O. (2004). RETRASO, a code for modeling reactive transport in saturated and unsaturated porous media. Geologica acta, 2(3), 235.
Samper, J., Juncosa, R., Delgado, J., & Montenegro, L. (2000). CORE2D: A code for nonisothermal water flow and reactive solute transport. Users manual version 2.: ENRESA Technical Publication 06/2000.
Samper, J., Xu, T., & Yang, C. (2009). A sequential partly iterative approach for multicomponent reactive transport with CORE2D. Computational Geosciences, 13(3), 301-316, doi:10.1007/s10596-008-9119-5.
Sanford, W. E., & Konikow, L. F. (1989a). Simulation of calcite dissolution and porosity changes in saltwater mixing zones in coastal aquifers. Water Resources Research, 25(4), 655-667, doi:10.1029/WR025i004p00655.
Sanford, W. E., & Konikow, L. F. (1989b). Simulation of calcite dissolution and porosity changes in saltwater mixing zones in coastal aquifers. Water Resources Research, 25(4), 655-667, doi:10.1029/WR025i004p00655.
Santoro, A. (2010). Microbial nitrogen cycling at the saltwater–freshwater interface. Hydrogeology Journal, 18(1), 187-202, doi:10.1007/s10040-009-0526-z.
Sasman, R. T. (1974). The Future of Ground-Water Resources in DuPage County. Ground Water, 12(5), 277-282, doi:10.1111/j.1745-6584.1974.tb03032.x.
Sauter, F. J., & Leijnse, A. (1993). METROPOL's User's Guide. National Institute of Public Health and Environmental Protection, Bilthoven, the Netherlands.
Schäfer, D., Schäfer, W., & Kinzelbach, W. (1998). Simulation of reactive processes related to biodegradation in aquifers: 1. Structure of the three-dimensional reactive transport model. Journal of Contaminant Hydrology, 31(1–2), 167-186, doi:http://dx.doi.org/10.1016/S0169-7722(97)00060-0.
Schäfer, W., & Therrien, R. (1995). Simulating transport and removal of xylene during remediation of a sandy aquifer. Journal of Contaminant Hydrology, 19(3), 205-236, doi:http://dx.doi.org/10.1016/0169-7722(95)00018-Q.
Scheidegger, A. E. (1974). The physics of flow through porous media: University of Toronto Press.
Schincariol, R. A., & Schwartz, F. W. (1990). An experimental investigation of variable density flow and mixing in homogeneous and heterogeneous media. Water Resources Research, 26(10), 2317-2329, doi:10.1029/WR026i010p02317.
Schmork, S., & Mercado, A. (1969). Upconing of Fresh Water—Sea Water Interface Below Pumping Wells, Field Study. Water Resources Research, 5(6), 1290-1311, doi:10.1029/WR005i006p01290.
Seetharam, S. C. (2003). An investigation o f the thermo/hydro/chemical/mechanical behavior of unsaturated soils. Ph.D., Cardiff University, Wales, UK,
Segol, G., & Pinder, G. F. (1976). Transient simulation of saltwater intrusion in southeastern Florida. Water Resources Research, 12(1), 65-70, doi:10.1029/WR012i001p00065.
Segol, G., Pinder, G. F., & Gray, W. G. (1975a). A Galerkin-finite element technique for calculating the transient position of the saltwater front. Water Resources Research, 11(2), 343-347, doi:10.1029/WR011i002p00343.
Segol, G., Pinder, G. F., & Gray, W. G. (1975b). A Galerkin-finite element technique for calculating the transient position of the saltwater front. Water Resources Research, 11(2), 343-347, doi:10.1029/WR011i002p00343.
Shao, H., Dmytrieva, S. V., Kolditz, O., Kulik, D. A., Pfingsten, W., & Kosakowski, G. (2009a). Modeling reactive transport in non-ideal aqueous–solid solution system. Applied Geochemistry, 24(7), 1287-1300, doi:http://dx.doi.org/10.1016/j.apgeochem.2009.04.001.
Shao, H., Kulik, D. A., Berner, U. R. S., Kosakowski, G., & Kolditz, O. (2009b). Modeling the competition between solid solution formation and cation exchange on the retardation of aqueous radium in an idealized bentonite column. Geochemical Journal, 43(6), e37-e42, doi:10.2343/geochemj.1.0069.
Shuster, E. T., & White, W. B. (1971). Seasonal fluctuations in the chemistry of lime-stone springs: A possible means for characterizing carbonate aquifers. Journal of Hydrology, 14(2), 93-128, doi:http://dx.doi.org/10.1016/0022-1694(71)90001-1.
Simmons, C. (2005). Variable density groundwater flow: From current challenges to future possibilities. Hydrogeology Journal, 13(1), 116-119, doi:10.1007/s10040-004-0408-3.
Simmons, C. T., Fenstemaker, T. R., & Sharp Jr, J. M. (2001). Variable-density groundwater flow and solute transport in heterogeneous porous media: approaches, resolutions and future challenges. Journal of Contaminant Hydrology, 52(1–4), 245-275, doi:http://dx.doi.org/10.1016/S0169-7722(01)00160-7.
Simmons, C. T., Narayan, K. A., & Wooding, R. A. (1999). On a test case for density-dependent groundwater flow and solute transport models: The Salt Lake Problem. Water Resources Research, 35(12), 3607-3620, doi:10.1029/1999WR900254.
Simpson, M. J., & Clement, T. P. (2003). Theoretical analysis of the worthiness of Henry and Elder problems as benchmarks of density-dependent groundwater flow models. Advances in Water Resources, 26(1), 17-31, doi:http://dx.doi.org/10.1016/S0309-1708(02)00085-4.
Šimůnek, J., & Suarez, D. L. (1994). Two-dimensional transport model for variably saturated porous media with major ion chemistry. Water Resources Research, 30(4), 1115-1133, doi:10.1029/93WR03347.
Singhal, B. B. S., & Gupta, R. P. (2010). Applied Hydrogeology of Fractured Rocks: Second Edition: Springer.
Singurindy, O., & Berkowitz, B. (2003a). Evolution of hydraulic conductivity by precipitation and dissolution in carbonate rock. Water Resources Research, 39(1), 1016, doi:10.1029/2001WR001055.
Singurindy, O., & Berkowitz, B. (2003b). Evolution of hydraulic conductivity by precipitation and dissolution in carbonate rock. Water Resources Research, 39(1), n/a-n/a, doi:10.1029/2001WR001055.
Singurindy, O., Berkowitz, B., & Lowell, R. P. (2004a). Carbonate dissolution and precipitation in coastal environments: Laboratory analysis and theoretical consideration. Water Resources Research, 40(4), W04401, doi:10.1029/2003WR002651.
Singurindy, O., Berkowitz, B., & Lowell, R. P. (2004b). Carbonate dissolution and precipitation in coastal environments: Laboratory analysis and theoretical consideration. Water Resources Research, 40(4), W04401, doi:10.1029/2003WR002651.
Souza, W. R., & Voss, C. I. (1987). Analysis of an anisotropic coastal aquifer system using variable-density flow and solute transport simulation. Journal of Hydrology, 92(1–2), 17-41, doi:http://dx.doi.org/10.1016/0022-1694(87)90087-4.
Spycher, N. F., Sonnenthal, E. L., & Apps, J. A. (2003). Fluid flow and reactive transport around potential nuclear waste emplacement tunnels at Yucca Mountain, Nevada. Journal of Contaminant Hydrology, 62–63(0), 653-673, doi:http://dx.doi.org/10.1016/S0169-7722(02)00183-3.
Steefel, C. I., & Lasaga, A. C. (1994). A coupled model for transport of multiple chemical species and kinetic precipitation/dissolution reactions with application to reactive flow in single phase hydrothermal systems. American Journal of Science, 294(5), 529-592, doi:10.2475/ajs.294.5.529.
Steefel, C. I., & Lichtner, P. C. (1998). Multicomponent reactive transport in discrete fractures: II: Infiltration of hyperalkaline groundwater at Maqarin, Jordan, a natural analogue site. Journal of Hydrology, 209(1–4), 200-224, doi:http://dx.doi.org/10.1016/S0022-1694(98)00173-5.
Steefel, C. I., & MacQuarrie, K. T. B. (1996). Approaches to modeling of reactive transport in porous media. Reviews in Mineralogy and Geochemistry, 34(1), 85-129.
Steefel, C. I., & Yabusaki, S. B. (1996). OS3D/GIMRT, Software for multicomponent-multidimensional reactive transport: User’s Manual and Programmer’s Guide. Pacific Northwest National Laboratory, Richland, Washington.
.
Stetzenbach, K. J., Farnham, I. M., Hodge, V. F., & Johannesson, K. H. (1999). Using multivariate statistical analysis of groundwater major cation and trace element concentrations to evaluate groundwater flow in a regional aquifer. Hydrological Processes, 13(17), 2655-2673, doi:10.1002/(SICI)1099-1085(19991215)13:17<2655::AID-HYP840>3.0.CO;2-4.
Stoessell, R. K., WARD, W. C., FORD, B. H., & SCHUFFERT, J. D. (1989). Water chemistry and CaCO3 dissolution in the saline part of an open-flow mixing zone, coastal Yucatan Peninsula, Mexico. Geological Society of America Bulletin, 101(2), 159-169, doi:10.1130/0016-7606(1989)101<0159:wcacdi>2.3.co;2.
Stumm, W., & Morgan, J. J. (1970). Aquatic chemistry: an introduction emphasizing chemical equilibria in natural waters (Vol. vol. 1): Wiley-Interscience.
Stumm, W., & Morgan, J. J. (1981). Aquatic chemistry: an introduction emphasizing chemical equilibria in natural waters: John Wiley.
Sun, Y., Petersen, J. N., & Clement, T. P. (1999). Analytical solutions for multiple species reactive transport in multiple dimensions. Journal of Contaminant Hydrology, 35(4), 429-440, doi:http://dx.doi.org/10.1016/S0169-7722(98)00105-3.
Taniguchi, M., Ishitobi, T., & Shimada, J. (2006). Dynamics of submarine groundwater discharge and freshwater-seawater interface. Journal of Geophysical Research: Oceans, 111(C1), C01008, doi:10.1029/2005JC002924.
Tartakovsky, A. M., Redden, G., Lichtner, P. C., Scheibe, T. D., & Meakin, P. (2008). Mixing‐induced precipitation: Experimental study and multiscale numerical analysis. Water Resources Research, 44(6).
Tebes-Stevens, C., J. Valocchi, A., VanBriesen, J. M., & Rittmann, B. E. (1998). Multicomponent transport with coupled geochemical and microbiological reactions: model description and example simulations. Journal of Hydrology, 209(1–4), 8-26, doi:http://dx.doi.org/10.1016/S0022-1694(98)00104-8.
Therrien, R., & Sudicky, E. A. (1996). Three-dimensional analysis of variably-saturated flow and solute transport in discretely-fractured porous media. Journal of Contaminant Hydrology, 23(1–2), 1-44, doi:http://dx.doi.org/10.1016/0169-7722(95)00088-7.
Thompson, C., Smith, L., & Maji, R. (2007). Hydrogeological modeling of submarine groundwater discharge on the continental shelf of Louisiana. Journal of Geophysical Research: Oceans, 112(C3), C03014, doi:10.1029/2006JC003557.
Thorne, D., Langevin, C. D., & Sukop, M. C. (2006). Addition of simultaneous heat and solute transport and variable fluid viscosity to SEAWAT. Computers & Geosciences, 32(10), 1758-1768, doi:http://dx.doi.org/10.1016/j.cageo.2006.04.005.
Tirado, M., Clarke, R., Jaykus, L., McQuatters-Gollop, A., & Frank, J. (2010). Climate change and food safety: A review. Food Research International, 43(7), 1745-1765.
Tirado, M. C., Clarke, R., Jaykus, L. A., McQuatters-Gollop, A., & Frank, J. M. (2010). Climate change and food safety: A review. Food Research International, 43(7), 1745-1765, doi:http://dx.doi.org/10.1016/j.foodres.2010.07.003.
Van Camp, M., Mtoni, Y., Mjemah, I. C., Bakundukize, C., & Walraevens, K. (2014). Investigating seawater intrusion due to groundwater pumping with schematic model simulations: The example of the Dar es Salaam coastal aquifer in Tanzania. Journal of African Earth Sciences, 96(0), 71-78, doi:http://dx.doi.org/10.1016/j.jafrearsci.2014.02.012.
Van Der Lee, J. (1998). Thermodynamic and mathematical concepts of CHESS. In U. G. S. Technical Report 99-4259, USA (Ed.): Ecole de Mines de Paris, Fontainbleau, France.
van der Lee, J., De Windt, L., Lagneau, V., & Goblet, P. (2003). Module-oriented modeling of reactive transport with HYTEC. Computers & Geosciences, 29(3), 265-275, doi:http://dx.doi.org/10.1016/S0098-3004(03)00004-9.
Vandenbohede, A., & Lebbe, L. (2006). Occurrence of salt water above fresh water in dynamic equilibrium in a coastal groundwater flow system near De Panne, Belgium. Hydrogeology Journal, 14(4), 462-472, doi:10.1007/s10040-005-0446-5.
Vandenbohede, A., & Lebbe, L. (2007). Effects of tides on a sloping shore: groundwater dynamics and propagation of the tidal wave. Hydrogeology Journal, 15(4), 645-658, doi:10.1007/s10040-006-0128-y.
Vandenbohede, A., & Lebbe, L. (2011). Heat transport in a coastal groundwater flow system near De Panne, Belgium. Hydrogeology Journal, 19(6), 1225-1238, doi:10.1007/s10040-011-0756-8.
Vandenbohede, A., Lebbe, L., Adams, R., Cosyns, E., Durinck, P., & Zwaenepoel, A. (2010). Hydrogeological study for improved nature restoration in dune ecosystems–Kleyne Vlakte case study, Belgium. Journal of Environmental Management, 91(11), 2385-2395, doi:http://dx.doi.org/10.1016/j.jenvman.2010.06.023.
Vandenbohede, A., Lebbe, L., Gysens, S., Delecluyse, K., & DeWolf, P. (2008a). Salt water infiltration in two artificial sea inlets in the Belgian dune area. Journal of Hydrology, 360(1–4), 77-86, doi:http://dx.doi.org/10.1016/j.jhydrol.2008.07.018.
Vandenbohede, A., Luyten, K., & Lebbe, L. (2008b). Effects of Global Change on Heterogeneous Coastal Aquifers: A Case Study in Belgium. Journal of Coastal Research, 160-170, doi:10.2112/05-0447.1.
Vandenbohede, A., Van Houtte, E., & Lebbe, L. (2008c). Groundwater flow in the vicinity of two artificial recharge ponds in the Belgian coastal dunes. Hydrogeology Journal, 16(8), 1669-1681, doi:10.1007/s10040-008-0326-x.
Volker, R. E., & Rushton, K. R. (1982). An assessment of the importance of some parameters for seawater intrusion in aquifers and a comparison of dispersive and sharp-interface modelling approaches. Journal of Hydrology, 56(3–4), 239-250, doi:http://dx.doi.org/10.1016/0022-1694(82)90015-4.
Voss, C., Simmons, C., & Robinson, N. (2010). Three-dimensional benchmark for variable-density flow and transport simulation: matching semi-analytic stability modes for steady unstable convection in an inclined porous box. Hydrogeology Journal, 18(1), 5-23, doi:10.1007/s10040-009-0556-6.
Voss, C. I. (1984a). A finite-element simulation model for saturated-unsaturated fluid-density-dependent ground-water flow with energy transport or chemically-reactive single-species solute transport. US Geol Surv Water Resour Invest 84-4369(3), 409, doi:10.1023/A:1006504326005.
Voss, C. I. (1984b). A finite-element simulation model for saturated-unsaturated fluid-density-dependent ground-water flow with energy transport or chemically-reactive single-species solute transport. US Geol Surv Water Resour Invest, 84-4369, 409pp.
Voss, C. I., & Souza, W. R. (1987). Variable density flow and solute transport simulation of regional aquifers containing a narrow freshwater-saltwater transition zone. Water Resources Research, 23(10), 1851-1866, doi:10.1029/WR023i010p01851.
Walter, A. L., Frind, E. O., Blowes, D. W., Ptacek, C. J., & Molson, J. W. (1994). Modeling of multicomponent reactive transport in groundwater: 1. Model development and evaluation. Water Resources Research, 30(11), 3137-3148, doi:10.1029/94WR00955.
Walther, M., Delfs, J. O., Grundmann, J., Kolditz, O., & Liedl, R. (2012). Saltwater intrusion modeling: Verification and application to an agricultural coastal arid region in Oman. Journal of Computational and Applied Mathematics, 236(18), 4798-4809, doi:http://dx.doi.org/10.1016/j.cam.2012.02.008.
Wang, Y., Guo, Q., Su, C., & Ma, T. (2006). Strontium isotope characterization and major ion geochemistry of karst water flow, Shentou, northern China. Journal of Hydrology, 328(3–4), 592-603, doi:http://dx.doi.org/10.1016/j.jhydrol.2006.01.006.
Watson, T. A., Werner, A. D., & Simmons, C. T. (2010). Transience of seawater intrusion in response to sea level rise. Water Resources Research, 46(12), W12533, doi:10.1029/2010WR009564.
Weatherill, D., Simmons, C. T., Voss, C. I., & Robinson, N. I. (2004). Testing density-dependent groundwater models: two-dimensional steady state unstable convection in infinite, finite and inclined porous layers. Advances in Water Resources, 27(5), 547-562, doi:http://dx.doi.org/10.1016/j.advwatres.2004.01.003.
Webb, M. D., & Howard, K. W. F. (2011). Modeling the Transient Response of Saline Intrusion to Rising Sea-Levels. Ground Water, 49(4), 560-569, doi:10.1111/j.1745-6584.2010.00758.x.
Wernberg, T. (1998). Multicomponent groundwater transport with chemical equilibrium and kinetics: model development and evaluation. Hydrological Sciences Journal, 43(2), 299-317, doi:10.1080/02626669809492123.
Werner, A. D., Bakker, M., Post, V. E. A., Vandenbohede, A., Lu, C., Ataie-Ashtiani, B., et al. (2013). Seawater intrusion processes, investigation and management: Recent advances and future challenges. Advances in Water Resources, 51(0), 3-26, doi:http://dx.doi.org/10.1016/j.advwatres.2012.03.004.
Werner, A. D., & Gallagher, M. R. (2006). Characterisation of sea-water intrusion in the Pioneer Valley, Australia using hydrochemistry and three-dimensional numerical modelling. Hydrogeology Journal, 14(8), 1452-1469, doi:10.1007/s10040-006-0059-7.
White, S. P. (1995). Multiphase Nonisothermal Transport of Systems of Reacting Chemicals. Water Resources Research, 31(7), 1761-1772, doi:10.1029/95WR00576.
White, W. B. (1969). Conceptual Models for Carbonate Aquifers. Ground Water, 7(3), 15-21, doi:10.1111/j.1745-6584.1969.tb01279.x.
Wicks, C. M., & Herman, J. S. (1996). Regional Hydrogeochemistry of a Modern Coastal Mixing Zone. Water Resources Research, 32(2), 401-407, doi:10.1029/95WR03244.
Wigley, T. M. L., & Plummer, L. N. (1976a). Mixing of carbonate waters. Geochimica et Cosmochimica Acta, 40(9), 989-995, doi:http://dx.doi.org/10.1016/0016-7037(76)90041-7.
Wigley, T. M. L., & Plummer, L. N. (1976b). Mixing of carbonate waters. Geochimica et Cosmochimica Acta, 40(9), 989-995, doi:http://dx.doi.org/10.1016/0016-7037(76)90041-7.
Wilson, A. M. (2005). Fresh and saline groundwater discharge to the ocean: A regional perspective. Water Resources Research, 41(2), W02016, doi:10.1029/2004WR003399.
Wolery, T. J. (1992). EQ3/6, A software package for geochemical modeling of aqueous systems: Package overview and installation guide. In U. G. S. Technical Report 99-4259, USA (Ed.): Lawrence Livermore National Laboratory, California.
Xu, T., & Li, J. (2013). Reactive Transport Modeling to Address the Issue of CO2 Geological Sequestration. Procedia Earth and Planetary Science, 7(0), 912-915, doi:http://dx.doi.org/10.1016/j.proeps.2013.03.153.
Xu, T., Samper, J., Ayora, C., Manzano, M., & Custodio, E. (1999). Modeling of non-isothermal multi-component reactive transport in field scale porous media flow systems. Journal of Hydrology, 214(1–4), 144-164, doi:http://dx.doi.org/10.1016/S0022-1694(98)00283-2.
Xu, T., Sonnenthal, E., Spycher, N., & Pruess, K. (2006a). TOUGHREACT-A simulation program for non-isothermal multiphase reactive geochemical transport in variably saturated geologic media: Applications to geothermal injectivity and CO2 geological sequestration. Comput. Geosci., 32(2), 145-165, doi:10.1016/j.cageo.2005.06.014.
Xu, T., Sonnenthal, E., Spycher, N., & Pruess, K. (2006b). TOUGHREACT—A simulation program for non-isothermal multiphase reactive geochemical transport in variably saturated geologic media: Applications to geothermal injectivity and CO2 geological sequestration. Computers & Geosciences, 32(2), 145-165, doi:http://dx.doi.org/10.1016/j.cageo.2005.06.014.
Yechieli, Y., Shalev, E., Wollman, S., Kiro, Y., & Kafri, U. (2010). Response of the Mediterranean and Dead Sea coastal aquifers to sea level variations. Water Resources Research, 46(12), W12550, doi:10.1029/2009WR008708.
Yechieli, Y., & Wood, W. W. (2002). Hydrogeologic processes in saline systems: playas, sabkhas, and saline lakes. Earth-Science Reviews, 58(3–4), 343-365, doi:http://dx.doi.org/10.1016/S0012-8252(02)00067-3.
Yeh, G. T., & Tripathi, V. S. (1989). A critical evaluation of recent developments in hydrogeochemical transport models of reactive multichemical components. Water Resources Research, 25(1), 93-108, doi:10.1029/WR025i001p00093.
Yeh, G. T., & Tripathi, V. S. (1990). HYDROGEOCHEM, A Coupled Model of HYDROlogic Transport and GEOCHEMical Equilibria in Reactive Multicomponent Systems.: ORNL-6371. Oak Ridge National Laboratory.
Yoon, H., Valocchi, A. J., Werth, C. J., & Dewers, T. (2012). Pore‐scale simulation of mixing‐induced calcium carbonate precipitation and dissolution in a microfluidic pore network. Water Resources Research, 48(2).
Younes, A., Ackerer, P., & Mose, R. (1999). Modeling Variable Density Flow and Solute Transport in Porous Medium: 2. Re‐Evaluation of the Salt Dome Flow Problem. Transport in Porous Media, 35(3), 375-394, doi:10.1023/A:1006504326005.
Younes, A., & Fahs, M. (2013). A semi-analytical solution for the reactive Henry saltwater intrusion problem. Water, Air, & Soil Pollution, 224(11), 1-11, doi:10.1007/s11270-013-1779-7.
Zhang, Q., Volker, R. E., & Lockington, D. A. (2004). Numerical investigation of seawater intrusion at Gooburrum, Bundaberg, Queensland, Australia. Hydrogeology Journal, 12(6), 674-687, doi:10.1007/s10040-004-0333-5.
Zheng, C., & Wang, P. P. (1999). A modular three-dimensional multi-species transport model for simulation of advection, dispersion and chemical reactions of contaminants in groundwater systems; documentation and user's guide. Contract report SERDP-99-1, U.S. Army Engineer Research and Development Center Vicksburg, MS.
Zimmermann, S., Bauer, P., Held, R., Kinzelbach, W., & Walther, J. H. (2006). Salt transport on islands in the Okavango Delta: Numerical investigations. Advances in Water Resources, 29(1), 11-29, doi:http://dx.doi.org/10.1016/j.advwatres.2005.04.013.
Zysset, A., Stauffer, F., & Dracos, T. (1994). Modeling of chemically reactive groundwater transport. Water Resources Research, 30(7), 2217-2228, doi:10.1029/94WR00230.
Copyright Notice
© Licențiada.org respectă drepturile de proprietate intelectuală și așteaptă ca toți utilizatorii să facă același lucru. Dacă consideri că un conținut de pe site încalcă drepturile tale de autor, te rugăm să trimiți o notificare DMCA.
Acest articol: جامعة تونس المنار [301438] (ID: 301438)
Dacă considerați că acest conținut vă încalcă drepturile de autor, vă rugăm să depuneți o cerere pe pagina noastră Copyright Takedown.