Impact of Aquifer Heterogeneity on Subsurface Flow and Salt Transport at Different Scales From a Method to Determine Parameters of Heterogeneous Permeability at Local Scale to a Large-Scale Model for the Sedimentary Basin of Thuringia Dissertation zur Erlangung des akademischen Grades doctor rerum naturalium (Dr. rer. nat.) vorgelegt dem Rat der Chemisch-Geowissenschaftlichen Fakultät der Friedrich-Schiller-Universität Jena von Dipl.-Math. Alraune Zech geboren am 02. Juni 1983 in Quedlinburg Gutachter: Prof.Dr. Sabine Attinger Friedrich-Schiller-Universität Jena Prof.Dr. Olaf Kolditz Technische Universität Dresden Tag der öffentlichen Verteidigung: 20. September 2013 Abstract Groundwater is an important natural resource. More than two thirds of the German population rely on groundwater as the daily drinking water supply. Hence, the understanding of processes like fluid flow and solute transport in porous media are of broad relevance. A profound process understanding is not only important for rational management of water resources but in particular for the preservation of subsurface water quality. Further relevant aspects are the optimization of irrigation and drainage efficiency, safe and economic extraction of subsurface mineral and energy resources, as well as subsurface storage of energy and waste. Subsurface hydrology faces two problems. On the one hand, information about the character of the subsurface structure and hydraulic properties is scarce. Data of quantities, like hydraulic pressure, porosity, and permeability, are limited to a few locations. It results in a lack of spatial resolution of known aquifer parameters for large-scale flow and transport problems. This is aggravated by the fact that the processes are very slow. Thus, there is a deficiency in the temporal resolution of observations. On the other hand, many aquifer parameters are heterogeneously distributed in space. In particular, permeability shows a spatial variability that results from the complex geological processes through which aquifers evolve. Large differences in the permeability can be observed on local scale, i.e. in the range of meters, as well as on large scale, i.e. in the range of kilometers up to the scale of geological structures as the Thuringian Basin. However, permeability is a very important hydraulic property since it controls the groundwater flow velocity and hence flow and salt transport. In general, the scarcity of data does not allow to give a detailed spatial resolution of heterogeneously distributed permeability. The coincidence of both issues, scarcity of information and heterogeneity, inhibits a simple analysis of the most problems of subsurface hydrology. This thesis is dedicated to the question of how aquifer heterogeneity impacts on processes of flow and salt transport in the subsurface at different scales. We present a large-scale numerical model of the Thuringian basin in order to investigate the mechanisms of brine transport within the aquifers of a shallow sedimentary basin. The emphasis lies on the effects that heterogeneous permeability have on the flow and salt patterns. To perform numerical simulations, prior knowledge of hydraulic parameters is necessary. Of particular interest are those parameters that describe the heterogeneous structure of porous media. Therefore, the second major issue of this work is the development and discussion of a method to determine the statistical parameters which describe the spatial distribution of permeability. The first part of the work refers to the basic principles of modeling flow and salt transport in porous media. We introduce the physical equations which describe the relevant subsurface processes of fluid flow, salt transport and heat transport. We explain the basic concepts of geostatistics and upscaling theory. The focus lies on the statistical description of permeability, which is the basis for the upscaling methods developed in the second part. In addition, we present numerical benchmarks in order to indicate the usability of the simulation software OpenGeoSys. We use this numerical solver for the modeling of both, the complex large-scale model of the sedimentary basin of Thuringia in the second part and for the pumping test model in the third part. In the second part, we present a numerical modeling approach to investigate the fluid dynamics in the Thuringian basin. We focus on the impact of aquifer heterogeneity and fluid density differences on brine transport. Central questions are: How does the large-scale fluid dynamics look like? Does a coupling between thermallyinduced deep fluid convection and near-surface groundwater flow exist? How comes that saline groundwater reaches or comes close to the surface, which is a phenomenon that is observed in many places in the Thuringian Basin? We carry out numerical simulations of fluid flow, mass and heat transport in order to understand the role of geological features such as faults, aquifer heterogeneity, as well as fluid density differences caused by temperature and salt concentration gradients. For this purpose we construct a profile model that represents the geological setting of the Thuringian basin incorporating major hydraulic units and fault structures. The numerical results indicate that regional groundwater flow determines brine migration. The pattern of groundwater flow depends strongly on the local hydraulic parameters. A qualitative sensitivity analysis indicates that small variations in permeability can have significant influence on the flow and salt patterns. The local mean value, the degree of heterogeneity, and the local correlation structure of permeability impact on the location and amount of dissolved salt. This directly affects the amount of salt which is transported to near-surface regions. Also variation in fluid density due to salt concentration differences can cause significant changes in the flow pattern. Brine is heavier than fresh water. If saline groundwater is transported upwards, e.g. by pressure induced groundwater flow, then the higher fluid density of brine causes a counteracting downward movement. This effect leads to enhanced mixing. Increased mixing amplifies the salinization of the deep aquifers but prevents the upward movement of highly concentrated brine. Contrariwise, the simulations show that temperature can be neglected as driving mechanism for fluid flow. The shallow basin structure inhibits the developments of thermal convection on a regional scale, due to small temperature differences. The third part of the thesis refers to an upscaling method which we develop in order to determine the statistics of aquifer heterogeneity from pumping test data. Permeability is the key parameter to describe groundwater flow velocities and the salt distribution. However, the effects of heterogeneity cannot be captured by a single mean value. According to the statistical approach, permeability can be described by a log-normal distribution which provides additional parameters, like the variance and the correlation length. The variance represents the degree of heterogeneity, whereas the correlation length describes the distance up to which permeability values are correlated in space. Head measurements of pumping tests are commonly used to estimate hydraulic properties of porous media. Therefore it is reasonable to develop analytical tools, which allow the determination of parameters of aquifer heterogeneity from pumping test data. First, we derive a formula for the hydraulic head describing the mean drawdown of a three dimensional steady state pumping test in heterogeneous anisotropic porous media. The derivation is based on the Coarse Graining upscaling method. The closed form solution of the effective well flow hydraulic head can be understood as an extension of Thiem’s formula to heterogeneous porous media. The effective well flow solution is a function of the radial distance and accounts for the statistics of the permeability, namely geometric mean variance, horizontal correlation length and anisotropy ratio. We exploit the nature of the analytical solution to perform a sensitivity analysis on the parameters of the effective well flow head solution and implement an inverse estimation strategy. We analyze numerical pumping tests, both an ensemble of as well as single pumping tests, to show the applicability of the head solution to interpret drawdown data. The results of the inverse estimation procedure show excellent agreement of estimated statistical parameters with initial values. Second, we determine whether the analytical solution of effective well flow is capable of providing accurate and confident parameter estimates of a heterogeneous aquifer under limited data availability. This is of practical relevance since head measurements are limited in on-site pumping tests. We use simulated pumping tests to systematically reduce sampling size while also determining the accuracy and uncertainty of estimates at each level of data availability. Our findings indicate that the accuracy and uncertainty of estimated parameters are sensitive to the number and spatial distribution of head measurements. Piezometers are required at both locations, i.e. directly at the pumping well and at large distance, to reliably estimate the respective values of hydraulic conductivity. Likewise, several piezometers are needed in the vicinity of the well to increase accuracy and reduce uncertainty in estimates of horizontal correlation length. We then apply the same analytical solution to estimate the statistical parameters of a fluvial heterogeneous aquifer at the test site Horkheimer Insel, Germany. The estimated mean conductivity, variance and correlation length agree very good with results from laboratory measurements. Our work provides valuable implications regarding the conceptual design of ground water pumping tests and the predictive power of established pumping test sites. 6 Zusammenfassung Grundwasser ist eine wichtige natürliche Ressource. In Deutschland wird der tägliche Bedarf an Trinkwasser von mehr als zwei Drittel der deutschen Bevölkerung durch Grundwasser gedeckt. Aus diesem Grund ist das Verständnis für die Prozesse im Untergrund wie Fluidfluss und Stofftransport von weitreichender Relevanz. Ein profundes Prozessverständnis ist sowohl für ein vernünftiges Ressourcenmanagement als auch für die Sicherung der Grundwasserqualität wichtig. Weitere relevante Aspekte sind die Optimierung von Bewässerung und Entwässerung, die sichere und ökonomische Entnahme von Bodenschätzen und Energie, sowie die Einlagerung von Energie und Abfällen im Untergrund. Die Hydrogeologie sieht sich zwei grundlegenden Problemen gegenüber. Zum einen gibt es oft sehr wenige Informationen zur Beschaffenheit des Untergrunds sowie zu hydraulischen Eigenschaften. Messdaten von physikalischen Größen, wie piezometrischer Druck, Porosität und Leitfähigkeit sind beschränkt auf wenige Orte. Dadurch ergibt sich ein Mangel an räumlicher Auflösung von bekannten Daten für großskalige Fließ-und Transportprobleme. Hinzukommt, dass die Prozesse im Untergrund sehr langsam ablaufen, was eine geringe zeitliche Auflösung der Beobachtungsdaten zur Folge hat. Das zweite Problem bezieht sich auf die räumliche Heterogenität vieler Aquiferparameter. Insbesondere die Leitfähigkeit unterliegt sehr großen Schwankungen sowohl auf kleiner Skala, d.h. im Bereich von einigen Metern, als auch auf großer Skala, d.h. im Bereich von Kilometern, bis hin zu Ausdehnungen von großräumigen geologischen Strukturen wie dem Thüringer Sedimentbecken. Die starke räumliche Heterogenität resultiert aus den komplexen geologischen Entwicklungsprozessen. Die Permeabilität wiederum ist eine sehr wichtige hydraulische Eigenschaft des Bodens, denn sie kontrolliert die Grundwasserfließgeschwindigkeit und damit die Fließmuster und die daraus resultierende Salzverteilung. Im Allgemeinen erlaubt der Mangel an Daten keine detailgetreue räumliche Auflösung der heterogenen Leitfähigkeiten. Das Vorhandensein beider Aspekte, Informationsmangel und Heterogenität, erschweren eine einfache Analyse der meisten hydrogeologischen Fragestellung maßgeblich. Die vorliegende Arbeit ist der Frage gewidmet wie Aquiferheterogenität Fließ-und Transportprozesse im Untergrund auf verschiedenen Skalen beeinflusst. Zum einen wird ein großräumiges numerisches Untergrundmodell des Thüringer Beckens verwendet um die Mechanismen des Salztransports in Aquiferen eines flachen Sedimentbeckens zu untersuchen. Der Schwerpunkt liegt auf den Effekten die heterogene Leitfähigkeiten auf die Fließ-und Salzmuster haben. Um numerischen Simulationen durchführen zu können, werden Werte für hydrologische Parameter benötigt. Von besonderem Interesse sind dabei solche Parameter, die die heterogene Struktur des Untergrunds beschreiben. Aus diesem Grund liegt der zweite Schwerpunkt dieser Arbeit auf der Entwicklung und Diskussion einer Skalierungsmethode. Diese dient der Bestimmung von statistischen Parametern, mit denen sich die räumliche Verteilung von Leitfähigkeiten beschreiben lässt, aus Pumptestdaten. Im ersten Teil der Arbeit werden die Grundlagen der Modellierung von Fließ- und Transportprozessen im Untergrund eingeführt. Beschrieben werden unter anderem die physikalischen Gleichungen der relevanten Prozesse Fluidfluss, Salz-und Wärmetransport. Weiterhin werden die grundlegenden Konzepte der Geostatistik sowie der Skalierungstheorie (Upscaling) erklärt. Der Fokus liegt auf der statistischen Beschreibung von heterogener Leitfähigkeit, welche Grundlage für die Methode ist, die im dritten Teil der Arbeit entwickelt wird. Zusätzlich werden numerische Benchmarks präsentiert um die Verwendbarkeit der Simulationssoftware OpenGeoSys zu zeigen. Sowohl für die Modellierung des Thüringer Beckens im zweiten als auch für die numerischen Pumptests im dritten Teil, wird dieser numerische Löser verwendet. Ein numerischer Modelansatz zur Untersuchung der Fluiddynamik des Thüringer Beckens wird im zweiten Teil der Arbeit präsentiert. Der Fokus liegt auf dem Einfluss von Aquiferheterogenitäten und der Fluiddichte auf den Salztransport. Zentrale Fragen sind: Wie sieht die großräumige tiefe und oberflächennahe Fluiddynamik aus? Existieren Kopplungen thermisch angetriebener tiefer Fluidkonvektion mit oberflächennahen Grundwasserströmungen? Wie ist das Phänomen zu erklären, dass salzige Grundwasser in oberflächennahen Schichten vorkommen, was an vielen Orten im Thüringer Becken zu beobachten ist? Numerische Simulationen zum Fluidfluss, Salztransport sowie zum Wärmetransport im Thüringer Becken werden durchgeführt um zu verstehen, welche Bedeutung geologische Besonderheiten wie Verwerfungen, Aquiferheterogenität und Dichteunter- schiede haben. Schwankungen in der Fluiddichte können dabei durch Temperatur- und Salzkonzentrationsgradienten entstehen. Zu diesem Zweck wird ein geologisches Profilmodell konstruiert, welches die geologischen Begebenheiten mit den Hauptgrundwasserleitern und -stauern sowie Verwerfungen im Thüringer Becken abbildet. Die Simulationsergebnisse zeigen, dass der regionale Grundwasserfluss die Haupttriebkraft für die Salzverteilung ist. Der Grundwasserfluss wiederum ist stark von den lokalen hydraulischen Eigenschaften abhängig. Eine qualitative Sensitivitätsanalyse zeigt, dass kleine Veränderungen in der Leitfähigkeit großen Einfluss auf die Fließmuster und die Salzverteilung haben können. Der lokale Mittelwert, die Stärke der Abweichungen vom Mittelwert sowie die räumliche Korrelationsstruktur der Leitfähigkeit beeinflussen die Verteilung sowie die Menge an gelöstem Salz im Grundwasser. Davon ist insbesondere die Menge an Salz abhängig, die in oberflächennahe Bereiche transportiert wird. Ebenso können salzkonzentrationsabhängige Variationen in der Fluiddichte die Fließmuster maßgeblich verändern. Salzige Grundwasser sind schwerer als frisches Grundwasser. Werden hochkonzentrierte Solen durch druckgetriebenen Grundwasserfluss aufwärts transportiert, dann bewirkt die höhere Fluiddichte eine entgegengesetzte Absinkbewegung die Durchmischung intensiviert. Die verstärkte Durchmischung wiederum erhöht die Versalzung der tiefen Grundwasserleiter, aber verhindert den Aufstieg von hochkonzentrierter Sole. Auf der anderen Seite zeigen die Simulationen auch, dass Temperatur als Triebkraft für konvektiven Fluidfluss vernachlässigt werden kann. Die flache Struktur des Beckens verhindert die Entwicklung von großräumiger thermischer Konvektion aufgrund der geringen Temperaturunterschiede. Der dritte Teil der Arbeit beschäftigt sich mit einer Skalierungsmethode, welche dazu dient die statistischen Parameter der heterogenen Leitfähigkeitsverteilung mittels Pumptestdaten zu bestimmen. Leitfähigkeit ist der Schlüsselparameter um Grundwasserfließgeschwindigkeiten und die Salzverteilung zu beschreiben. Allerdings benötigt es mehr als einen Mittelwert um die Effekte von Heterogenitäten zu beschreiben. Bei Verwendung eines statistischen Ansatzes lassen sich Leitfähigkeiten durch eine Log-Normal-Verteilung beschreiben. Zusätzliche statistische Parameter wie die Varianz und die Korrelationslänge erlauben eine detailliertere Charakterisierung. Dabei bestimmt die Varianz die Stärke der lokalen Schwankungen, d.h. den Grad der Heterogenität. Die Korrelationslänge beschreibt die Distanz über die Leitfähigkeitswerte räumlich korreliert sind. Pumptests werden häufig verwendet um die hydraulischen Eigenschaften eines Aquifers zu bestimmen. Deshalb ist es sinnvoll analytische Methoden zu entwickeln, welche es erlauben die statistischen Parameter der heterogenen Leitfähigkeitsverteilung aus Pumptestdaten zu bestimmen. Zunächst wird eine Formel für das hydraulische Druckpotential hergeleitet, die den mittleren Absenktrichter eines kleinskaligen stationären Pumptests in einem dreidimensionalen, anisotropen, heterogenen Medium beschreibt. Die Herleitung basiert auf der Skalierungsmethode Coarse Graining. Die analytische Lösung für den effektiven Brunnenfluss kann man verstehen als eine Erweiterung für Thiems Formel auf heterogene Medien. Die Lösung des effektiven Brunnenflusses ist einerseits eine Funktion vom Abstand zum Pumpbrunnen und anderseits abhängig von den statistischen Größen Mittelwert, Varianz, Korrelationslänge und Anisotropierate der Leitfähigkeit. Die geschlossen-analytische Form der Lösung erlaubt es eine Sensitivitätsanalyse bezüglich der statistischen Parameter durchzuführen sowie eine inverse Schätzmethode zu entwickeln. Die Analyse von numerischen Pumptests, sowohl einzelner als auch von Ensembles, zeigt, dass Absenktrichter aus Pumptests in heterogenen Medien sehr genau durch die Lösung des effektiven Brunnenflusses beschrieben werden können. Die Ergebnisse der Inversen Schätzung der statistischen Leitfähigkeitsparameter stimmen sehr gut mit den Anfangswerten überein. Es folgt die Untersuchung, ob die analytische Lösung des effektiven Brunnenflusses dazu in der Lage ist genaue und zuverlässige Schätzwerte für die statistischen Parameter unter reduzierter Datenverfügbarkeit zu liefern. Dies ist relevant, da in der Praxis meist nur wenige Messungen des hydraulischen Druckpotentials verfügbar sind. Numerische Pumptest werden generiert und ausgewertet mit Hinblick auf die Zuverlässigkeit und die Genauigkeit der Schätzergebnisse unter reduzierter Datenverfügbarkeit. Dazu wird die Anzahl von Messdatenpunkten systematisch reduziert und für jedes Level der Datenverfügbarkeit die statistischen Parameter und ihre Unsicherheiten invers geschätzt. Die Ergebnisse zeigen, dass die Genauigkeit und Zuverlässigkeit der geschätzten statistischen Parameter sehr sensitiv gegenüber Anzahl und räumlicher Verteilung der Druckpotentialmessungen sind. Untersuchungsbrunnen sind sowohl im Nahfeld, d.h. in direkter Umgebung des Pumpbrunnens, als auch im Fernfeld notwendig um eine sinnvolle Schätzung von Heterogenitätsparametern gewährleisten zu können. Insbesondere die Schätzbarkeit der Korrelationslänge erfordert eine gewisse Anzahl an Piezometern im Umfeld des Pumpbrunnens. Abschließend wird die Schätzmethode auf Messdaten von Pumptests auf der Horkheimer Insel, Deutschland angewendet. Es wird gezeigt, dass die Schätzwerte der Methode des effektiven Brunnenflusses für die mittlere Leitfähigkeit, Varianz und Korrelationlänge sehr gut mit den Messergebnissen aus Laboruntersuchungen übereinstimmen. Die Untersuchungen geben wertvolle Hinweise in Bezug auf das konzeptionelle Design von Pumptests als auf die Güte der Schätzung von Parametern aus vorliegenden Pumptestdaten. List of Symbols and Abbreviations Symbol Unit Quantity l [m] longitudinal solute dispersivity t [m] transversal solute dispersivity ¯l [m] longitudinal thermal dispersivity ¯t [m] transversal thermal dispersivity C [-] solute expansion coefficient T [K-1] thermal expansion coefficient (e) anisotropy function G exponential integral function d Kronecker delta . [-] proportionality factor . correlation function ˜Fourier transform of correlation function . [mD] permeability . [m] filter constant s [kgm2 s-3 K-1] thermal conductivity of solid w [kgm2 s-3 K-1] thermal conductivity of water . [-] weighting factor µ mean of log-conductivity . [kgm-1 s-1] viscosity f [-] porosity . [-] angular coordinate . [kgm-3] fluid density 0 [kgm-3] fresh water density s [kgm-3] density of solid rock 2 [-] variance of log-conductivity ˆ2 [-] estimated variance of log-conductivity t [-] tortuosity . abbreviation for ln(Kefu/Kwell) w [m2] surface of the well V nabla operator A [-] relative space fraction of an element in a numerical mesh c [m2 s-2 K-1] thermal capacity of water cs [m2 s-2 K-1] thermal capacity of solid rock C [-] concentration C0, C1, C2 integration constants ¯ C [-] spatial mean concentration C [%] relative concentration difference Cr [-] CV d [-] D [m2 s-1] Dm [m2 s-1] Ds [m2 s-1] DT [m2 s-1] e [-] eˆ[-] f f˜ g/ g [ms-2] h [m] hˆ [m] hefw(r) [m] hh(r)) [m] hh(r)) [m] hThiem(r) [m] H [-] Kf [ms-1] KA [ms-1] KG [ms-1] KH [ms-1] Keq [ms-1] Kef [ms-1] Kefu [ms-1] Kwell [ms-1] K(x) [ms-1] KCG(r) [ms-1] Kˆ A [ms-1] Kˆ G [ms-1] Kˆ H [ms-1] Kˆ well [ms-1] Kˆ efu [ms-1] hKA) [ms-1] hKH) [ms-1] hKwell) [ms-1] i [m] `h [m] `v [m] `ˆ [m] L/Lz [m] Courant number covariance function number of spatial dimensions hydrodynamical dispersion tensor molecular diffusion coefficient mechanical dispersion tensor thermal diffusion/dispersion tensor anisotropy ratio estimated anisotropy ratio filter function Fourier transform of filter function gravitational acceleration / vector hydraulic head estimated hydraulic head effective well flow hydraulic head simulated hydraulic head ensemble mean of simulated hydraulic heads hydraulic head of Thiem’s solution Hurst coefficient hydraulic conductivity arithmetic mean conductivity geometric mean conductivity harmonic mean conductivity equivalent conductivity effective conductivity effective conductivity for uniform flow effective conductivity at the well heterogeneous conductivity distribution radial depending mean Coarse Graining conductivity estimated arithmetic mean conductivity estimated geometric mean conductivity estimated harmonic mean conductivity estimated conductivity at the well estimated conductivity far from field ’measured’ arithmetic mean conductivity ’measured’ geometric mean conductivity local mean of conductivity values at the well correlation length horizontal correlation length vertical correlation length estimated correlation length aquifer thickness/vertical domain extend ODE p PD PDE pdf Peg q q Q Qw r rw R Rac RaS RaT REF s S SRF SV t T ¯ T T Tf Z(x) [kgm-1 s-2] [m] [-] [ms-1] [ms-1] [m3 s-1] [m] [m] [m] [-] [-] [-] [m] [m3 kg-1] [s] [C] or [K] [-] [%] [m2 s-1] [ms-1] ordinary differential equation pressure penetration depth partial differential equation probability density function grid Peclet number Darcy velocity vector absolute Darcy velocity sink/source term pumping rate at the well radial distance radius of the well fixed radial distance from the pumping well critical Rayleigh number solute Rayleigh number thermal Rayleigh number representative elementary volume two-point separation distance specific storage coefficient spatial random function semivariogram function time temperature spatial mean temperature relative temperature difference transmissivity pore velocity vector spatial random function 14 Contents Abstract 3 Zusammenfassung 7 List of Symbols and Abbreviations 11 I. Modeling Fluid Flow and Salt Transport in a Sedimentary Basin 17 1. Principles of Fluid Flow and Salt Transport in Heterogeneous Porous Media 19 1.1. Introduction...................................... 19 1.2. PhysicsofPorousMedia............................... 20 1.3. GeostatisticsinSubsurfaceHydrology ....................... 25 1.4. Upscaling ....................................... 29 2. Numerical Benchmarks 35 2.1. Introduction...................................... 35 2.2. PumpingTestinHomogeneousPorousMedia . . . . . . . . . . . . . . . . . . . 35 2.3. UniformFlowinHeterogeneousPorousMedia. . . . . . . . . . . . . . . . . . . 37 2.4. TheElderBenchmark ................................ 40 2.5. TheSaltdomeProblem................................ 44 II. Impact of Heterogeneity and Density Dependent Flow on Regional Scale 51 3. Mechanisms of Salt Transport in the Thuringian Basin 53 3.1. Introduction...................................... 53 3.2. Hydrogeological Characterization of the Study Area . . . . . . . . . . . . . . . 55 3.3. NumericalGroundwaterFlowModel ........................ 59 3.4. DiscussionofSimulationResults .......................... 64 3.5. Conclusions ...................................... 84 III. Determining Aquifer Heterogeneity on Local Scale 87 4. The Extended Thiem’s Solution -Including the Impact of Heterogeneity 89 4.1. Introduction...................................... 89 4.2. StatementoftheProblem .............................. 91 4.3. MethodofCoarseGraining ............................. 93 4.4. TheExtendedThiem’sSolution........................... 95 4.5. InterpretingNumericalPumpingTest ....................... 99 4.6. SummaryandConclusions.............................. 105 5. Estimating Parameters of Aquifer Heterogeneity Using Pumping Tests 109 5.1. Introduction...................................... 109 5.2. TheoreticalFramework................................ 110 5.3. Pumping Tests under Limited Data Availability -a Paradigm for Field Applications113 5.4. FieldApplication................................... 118 5.5. Conclusions ...................................... 121 6. Summary and Outlook 123 A. Appendix 125 A.1. The Coarse Graining Conductivity in Anisotropic Media . . . . . . . . . . . . . 125 A.2.TheAnisotropyFunction .............................. 126 A.3. Derivation of the Effective Well Flow Head in Three Dimensions . . . . . . . . 127 A.4. Derivation of the Effective Well Flow Head in Two Dimensions . . . . . . . . . 128 Bibliography 131 Part I. Modeling Fluid Flow and Salt Transport in a Sedimentary Basin 18 1. Principles of Fluid Flow and Salt Transport in Heterogeneous Porous Media 1.1. Introduction During the last 20 years numerical modeling has gained significant importance in hydrogeology. The rapid evolution of computer technology allowed the development of complex software which reproduces the physical processes occurring in the subsurface with spatial and temporal resolution. Numerical models have several advantages. They can enhance the understanding of fluid dynamics due to the fact that they allow simultaneous evaluation of processes and thus investigation of coupling effects. Furthermore, models can serve quantitative aspects by providing estimates of quantities that can be used for engineering and management purposes. It requires several steps to establish a complex numerical model for subsurface flow and transport. Starting point is the formulation of a hydrogeological conceptual model. This aspect is significantly governed by the objective of the modeling project. The step contains the selection of the geological setup by identifying boundaries of the formation and significant geological features. The following step is to grasp the main mechanisms involved and to formulate physical equations describing the processes. In case existing numerical solvers are used, the third step contains the adaption of the hydrogeological conceptual model to the numerical framework. It requires the discretization of the domain in space and in time, to fix boundary conditions, initial conditions, as well as to set physical parameters. Developing a conceptual model is very specific to the problem under consideration. However, the physical equations which describe the major processes occurring in the subsurface are well known. We discuss them briefly in section 1.2. These conservation laws form the basis of physically-based simulation software, like OpenGeoSys [Kolditz et al., 2012b], FEFLOW [Fef, 2012], MODFLOW [U.S. Geological Survey, 2012], HydroGeoSphere [Brunner and Simmons, 2012], SUTRA [Voss and Provost, 2010], and many more. Complex numerical models with spatial and temporal resolution for thermal, hydraulic, and mechanical subsurface processes are important for the analysis of deep geological systems under high temperature, pressure and stress conditions. Application areas are e.g. reproduction and prediction of the solute plume development, geothermal energy utilization, nuclear waste disposal, and carbon dioxide storage in the deep geological formation. The hydraulic parameter of major interest within this work is the permeability. It shows a spatially heterogeneous distribution, which inhibits a simple analysis of many subsurface flow and transport problems. A detailed spatial resolution of the heterogeneous distribution is not possible in numerical models in general due to complexity. Therefore methods have been developed to handle aquifer heterogeneity in a stochastic manner which allows to incorporate uncertainty about the spatial distribution of permeability. We discuss the fundamental concepts of stochastic subsurface hydrology in section 1.3. To capture the aspects of heterogeneous permeability in large scale hydrogeological models, it is necessary to reduce complexity by a certain averaging strategy. This procedure is in general denoted upscaling. Upscaling methods 1. Principles of Fluid Flow and Salt Transport in Heterogeneous Porous Media aim to find representative parameters for spatially distributed quantities, which are capable of reproducing some average behavior of the system. We discuss the basics of upscaling in section 1.4. 1.2. Physics of Porous Media The simulation software of choice in this thesis is OpenGeoSys [Kolditz et al., 2012a, b]. It follows the basic idea of continuum mechanics that the evolution of a physical system is completely determined by conservation laws. Basic properties such as mass, momentum and energy are conserved during the considered process at all times. Any physical system can be completely determined by these conservation properties. To reproduce processes in porous media numerically, additional information concerning the consistencies of the material (e.g. fluids, solids, porous medium) in the form of constitutive laws is necessary. In the following we shortly introduce the continuum concept of porous media and state the macroscopic balance equations of mass of fluid, mass of solute, fluid momentum and energy conservation. It is not aimed to give derivations nor to discuss or expand the fundamentals of porous media theory. A detailed description of the theoretical background with derivations can be found for instance in Bear [1972], Kinzelbach and Rausch [1995], as well as Nield and Bejan [1999]. Detailed discussion on density dependent flow effects can be found e.g. in Kolditz et al. [1998] and Diersch and Kolditz [2002]. 1.2.1. Porous Media as Continuum Porous media is the geological material in the subsurface which consist of a solid matrix and pore space filled by water. The microscopic distribution of solid and pore space is complex and in general impossible to resolve. A description of the processes of flow and transport on microscopic level is neither feasible nor useful in most problems on macroscopic scale. The continuum approach targets at resolving these difficulties by a transition to the macroscopic level. At the macroscopic scale, porous media is regarded as a continuum that consists of an immobile solid phase and a mobile fluid phase. Both phases are characterized by macroscopic measurable quantities, which result from microscopic quantities by averaging over volumes. The transformation from microscopic to macroscopic level requires appropriate averaging strategies. The standard approach of averaging is to sample over a representative elementary volume (REV) [Bear, 1972]. On the one hand, the REV has to be sufficiently large for fluctuations of spatially averaged properties to be negligible on microscopic scale. On the other hand, it has to be small enough to be regarded as a point at the macroscopic scale. In hydrogeology the REV is the smallest volume over which measurements can be made that yield a representative value for the entire domain under consideration. Quantities and parameters of porous media are only valid above the scale of the REV. 1.2. Physics of Porous Media 1.2.2. Continuity Equation of Flow The continuity equation is derived from the principle of mass conservation of the fluid. It describes the evolution of fluid flow in space and time. It relates the temporal change in fluid mass to the differences in spatial distribution in combination with the existence of a fluid sink or source. The macroscopic mass balance equation of the fluid, averaged over a representative elementary volume (REV), in a porous medium is given by . (S) + r· (q)= Q. , (1.1) @t where S is the specific storage coefficient, f is the porosity, . is the fluid density, t is the time, q is the Darcy velocity vector, and Q. is the sink/source term of the fluid mass in an aquifer. 1.2.3. Darcy’s Law Darcy’s law characterizes the velocity vector q as a combination of porous medium characteristics and the hydraulic gradient. Darcy’s law can be understood as a momentum balance equation in porous media. In terms of hydraulic pressure p it is given by q = v = - . (rp + g) , (1.2) . where v is the fluid or pore velocity vector, . is the permeability of a porous medium, . is the viscosity, and g is the gravity vector. An alternative way to write Eq.(1.2) is the formulation of Darcy’s law with the freshwater hydraulic head h that directly relates the flow velocity to the driving forces q = -Kf (rh +(. - 0)/0ez) , (1.3) where Kf is the hydraulic conductivity, 0 is the fresh water density, . is the actual fluid density and ez is the unit vector in the gravitational direction. The buoyancy term (. - 0)/0ez is only present in processes where density differences are taken into account. In constant density processes the velocity is determined by the gradient in hydraulic head rh. The hydraulic conductivity 0g Kf = (1.4) . is a porous medium property since it depends on permeability . as a property of the solid phase and on the fluid characteristics density 0 and viscosity . The relation between hydrostatic pressure p and hydraulic head h is given by p h =+ zo, 0g where z0 denotes the vertical distance to the zero pressure position. 1. Principles of Fluid Flow and Salt Transport in Heterogeneous Porous Media 1.2.4. Solute Transport Equation The transport of a solute like salt in groundwater is governed by the advection-dispersion equation . (C) + r· (qC) -r· (D ·rC)= QC , (1.5) @t where C . [0, 1] is the relative mass concentration, D hydrodynamic dispersion tensor, and QC is the mass supply. Several fluxes are involved in the mass transport process in Eq. (1.5): the first term describes the temporal change, the second term describes an advective flux and the third term includes diffusive and dispersive effects. The right hand side accounts for a source of solute mass. The dispersion tensor D in Eq.(1.5) incorporates dispersive and diffusive fluxes. According to Scheidegger’s law [Bear, 1972] it is given in general tensor form by )qiqj Dij = Dm + Ds = Dm + tqij +( l - t, (1.6) q where t is the tortuosity and Dm is the coefficient of molecular diffusion. Ds is the mechanical dispersion with the longitudinal and transversal dispersivities l and t, ij is the Kronecker delta, q = |q| is the absolute value of Darcy velocity, and qi and qj are the Darcy velocities in ith and jth principal directions, respectively. A simplified description can be given in a Cartesian coordinate system, which is suitable for instance for uniform flow. Choosing the x1-axis to be parallel to the average Darcy velocity q the tensor Ds simplifies to 1 . BB . lq 00 0 tq 0 CCA Ds = . (1.7) 00 tq Thus the dispersion in main flow direction is given by the product of the average Darcy velocity q with the longitudinal dispersivity l. In the perpendicular directions the dispersion is proportional to the transversal dispersivity t, which is commonly assumed to be at least one order of magnitude smaller than l. According to Bear [1972], the axes of the coordinate system in which Ds is expressed in Eq. (1.7) are called the principal axes of dispersion. 1.2.5. Heat Transport Equation Energy conservation has to be taken into account when describing temperature dependent processes. The heat balance equation in porous media is based on a thermal equilibrium between solid and liquid phase, which can be assumed due to the low velocities. The equation is given by ss)@T (c. + (1 - )c + cq ·rT -r· (cDT ·rT )= QT , (1.8) @t where c. + (1 - )css is the heat capacity of the porous medium, with porosity , density . and specific heat capacity c of water as well as density s and specific heat capacity cs of 1.2. Physics of Porous Media solid. T is the temperature, t it the time, q is the Darcy velocity, and DT is the thermal diffusion/dispersion tensor. The heat source term is given in detail by QT = Qw + (1 - T )sQs , where Qw is the heat production of water and Qs is the heat production of solid. TT T The thermal diffusion-dispersion tensor DT is given in analogy to the hydrodynamic dispersion tensor D in Eq.(1.6), DT )qiqj ij = DT ij +¯ tqij +(¯ l - ¯t, (1.9) q w+(1-)s where DT = is the thermal diffusivity with heat conductivity w and s of water c. and solid, respectively. The second term and third term are part of the thermal dispersion tensor with heat dispersion coefficients ¯t and ¯l in transversal and longitudinal directions. ij is the Kronecker-delta, q = |q| is the absolute value of Darcy velocity, qi and qj denote the Darcy velocities in ith and jth principle directions, respectively. Due to the low velocities q in porous medium, thermal diffusion dominates the process of thermal dispersion. Therefore the thermal dispersion term in Eq. (1.9) is mostly neglected, resulting in cDT = w + (1 - )s in Eq.(1.8). 1.2.6. Equation of State for Fluid Density The description of processes where changes in fluid density take place requires a constitutive fluid density model. Based on experiments, equations are formulated which describe the relationship between fluid density and other relevant quantities like water pressure, solute concentration and temperature. Different descriptions exist, like an exponential function as given by Kolditz et al. [1998] or the additive description of Herbert et al. [1988]. The most common approximate relationship is the linearized equation of the bulk fluid density [Kolditz et al., 1998; Diersch and Kolditz, 2002; Nield and Bejan, 1999; Magri, 2005] accounting for thermal and solute mass effects . = 0 (1 + C C - T T ) , (1.10) where . denotes the fluid density, 0 is the reference density of freshwater, C . [0, 1] is the T -T0 change in relative solute concentration, C is the solute expansion coefficient, T = 2 T0 [0, 1] is the change in temperature relative to the reference temperature T0, and T is the thermal expansion coefficient. A similar model could be assumed for fluid viscosity since this property also depends on temperature and solute mass. However, investigations of Magri [2005] and Kaiser et al. [2011] show that the effects of fluid viscosity on flow pattern at basin scale is much less pronounced in comparison to density effects. Thus a constant viscosity model is appropriate for moderate temperature and solute mass gradients. 1. Principles of Fluid Flow and Salt Transport in Heterogeneous Porous Media 1.2.7. Stability Criteria for Density Dependent Flow The Rayleigh number is a classical instrument to characterize density dependent flow regimes. The value of the Rayleigh number characterizes the stability of the system and is used to predict the onset of convection. It results from a dimensional analysis of Darcy’s law (1.3) in combination with the equation of state for fluid density (1.10) and the transport equations for solutes (1.5) or heat (1.8). The Rayleigh number is defined with reference to the quantity of influence on the fluid density as Rayleigh number of solute RaS and thermal Rayleigh number RaT by Kf Lz C CKf Lz| T |T RaS = and RaT = , (1.11) D DT where Kf denotes the hydraulic conductivity, Lz is the vertical domain extent, f is the porosity, C is change in solute concentration, C is the solute expansion coefficient, D is the hydrodynamic dispersion in the direction of interest, T is the change in temperature, T is the thermal expansion coefficient, and DT is thermal diffusion-dispersion. In double diffusive phenomena, where both concentration and temperature differences cause buoyancy-driven transport, the Rayleigh number results as sum or difference of RaS and RaT, corresponding to the direction of impact [Nield and Bejan, 1999]. Perturbation theory identified two important critical Rayleigh numbers: Ra(1) =4p which characterizes the onset of c convection in stable stationary rolls and Ra(2) ˜ 400 marking the limit between stable and c unstable patterns. 1.2.8. Quantitative Indicators Characteristic measures can help to quantify differences in simulation results. This issue is in particular important for benchmarking and for sensitivity analysis of parameters. Indicators can underline visually observed results or quantify discrepancies between simulations which cannot be detected by pure visual comparison. A quantitative indicator has to be chosen according to the problem of interest. For instance Prasad and Simmons [2005] presented a number of methods to compare simulation results for variable density flow problems. We distinguish between indicators referring to a single simulation and comparative measures which quantify the difference between two simulations. Typical representatives for the first group are spatial and/or temporal means as integrated measures. Further characteristics are the penetration depth of concentration isolines, the center of mass, minimal or maximal velocities, breakthrough curves, and so forth. Indicatiors of the second group compare results which requires the analysis of the spatial structure and the distribution of quantities for of different simulations. This can be a quite challenging task. Various disciplines use pattern matching algorithms, like image analysis, meteorology or geographical information science. For a broad literature scan see e.g. Hagen-Zanker [2006]. We focus on two integrated measures. 1.3. Geostatistics in Subsurface Hydrology Mean Concentration We define the mean for the concentration C as the spatial average over the entire domain of a simulation k with an arbitrary shape of the numerical mesh by N X ¯ Ck = A(i)Ck(i) and (1.12) i=1 where N is the number of elements of the mesh, A(i) . [0;1] is the relative space fraction of ¯ the element i withPNi=1 A(i) = 1; C(i) is the concentration value of the element i. C is a dimensionless quantity due to the fact that the concentration is given in this work in relative ¯ amounts C . [0;1]. The mean concentration C corresponds directly to the amount of salt being solved within the domain. Relative Concentration Difference and Relative Temperature Difference To perform a comparison between two simulations k and l on the same numerical grid, we define the relative difference in concentration and temperature according to N XP Ckl = 100% A(i) |Ck(i) - Cl(i)| and i=1 (1.13) N 100% XP T kl = A(i) |Tk(i) - Tl(i)| , Tmax - Tmin i=1 where N is the number of elements of the mesh, A(i) . [0;1] is the relative space fraction of the element i withPNi=1 A(i) = 1. Ck(i) and Cl(i) are the concentration values and Tk(i) and Tl(i) are the temperature of element i of the simulations k and l, respectively. The division by the difference between the minimal temperature Tmin and maximal temperature Tmax is necessary to gain a relative quantity T . The relative difference measures the accordance of spatial pattern between two simulations. A value of zero corresponds to similarity. 1.3. Geostatistics in Subsurface Hydrology Subsurface hydrology faces two problems concerning the characterization of hydraulic properties of porous media: a strong variability of properties in space and the scarcity of information. The coincidence of both issues complicates the description and prediction of processes like the fluid, mass and heat transport in porous media. To account for the heterogeneity of subsurface characteristics, the theory of geostatistics has been developed in the context of hydrogeology in the 1960s. One of the first scientists using this concept was Matheron [1967]. Geostatistics attempts to find statistical characterizations which capture the patterns of spatial variability and to explain field observations. It further aims to recognize and quantify the impact of subsurface heterogeneity on flow and transport processes. Variables of interest are the hydraulic conductivity, transmissivity, hydraulic head, concentration, fracture density, and dispersivity. 1. Principles of Fluid Flow and Salt Transport in Heterogeneous Porous Media A detailed discussion on multiple geostatistical aspects are given e.g. in Dagan [1989]; Gelhar [1993] and Rubin [2003]. In the following we focus on hydraulic conductivity Kf . All statements and results are equally valid for the porous media quantity permeability, which is in one to one correspondence to the hydraulic conductivity trough the relation given in Eq. (1.4). The spatial variability of porous media is a result of complex geological processes through which aquifers evolve. Physical and chemical processes, like structural deformation, deposition and diagenesis control the geometry and texture of sedimentary deposits [Miall, 1996]. Field studies indicate that hydraulic conductivity can vary by orders of magnitude over short distances of the order of meters [Sudicky, 1986; Gelhar, 1993]. The spatial variability of hydraulic properties influences fluid flow and transport, due to the fact that the conductivity controls direction and magnitude of flow by Darcy’s law. Furthermore, heterogeneity triggers mechanical dispersion and dominates the movement of a solute plume. However, it is of practical significance to understand the actual movement of contaminants for their removal by remediation technologies. A deterministic approach to characterize complex spatial patterns is inappropriate because it requires the estimation of a large number of parameters which is in contrast to the rather small number of measurements available. The lack of perfect or complete measurements makes our knowledge of the variable uncertain , thus justifying the probabilistic approach. 1.3.1. Concept of Spatial Random Function The geostatistical description of a spatially distributed quantity, like hydraulic conductivity, uses the concept of spatial random functions (SRF). These are random variables that depend on the location and exhibit a stochastic spatial structure. The goal of constructing a subsurface quantity as SRF is to reduce the ensemble of measurements to a few useful statistics which capture mathematically the pattern of spatial variability. A SRF Z(x) can be characterized by a certain set of statistical parameters, in general its statistical moments. The spatial variation of a SRF often shows an overall trend and a random component associated with erratic fluctuations. These characteristics can be captured by the first two moments of Z(x). The first order moment is the expected value or arithmetic mean ZA(x), marking the trend. The second order moments include variance Var[Z(x)], the auto-covariance CV[Z(x)] and the semivariogram SV[Z(x)]. The interdependency of conductivity values in space is covered by the model of spatial correlation, which is expressed by the covariance function or semivariogram. In linear geostatistics higher moments are neglected. For hydraulic conductivity the property of stationarity is often assumed, meaning that the statistics do not change over space. Thus the mean is constant and the auto-correlation does not depend on x, but on the separation distance of two points s. Additionally, there exists a one-on-one correspondence between auto-covariance CV and semivariogram SV: ZA = E[Z(x)], CV(s) = Cov[Z(x + s),Z(x)] = E[Z(x + s)Z(x)] - ZA2 , SV(s) = Var[Z(x)] - CV(s). 1.3. Geostatistics in Subsurface Hydrology 1.3.2. Hydraulic Conductivity as Log-normally Distributed SFR It is common to model the hydraulic conductivity K(x) conceptually as log-normally distributed SRF, referring to the summary of several field studies given by Gelhar [1993]. It means that Y (x) = ln K x) is normally distributed with a Gaussian probability density func ( (  1(x-µ)2 tion pdfY (x)= 22 exp -in uni-variate form with µ and 2 being the mean and the 22 variance of Y , respectively. The probability density function of K(x) is given by  1 ln(x) - µ pdfK(x)= v exp-. (1.14) 22 x 22 The moments of K(x) can be calculated using the statistical parameters µ and 2. The arithmetic mean KA (first moment), the geometric mean KG as well as the harmonic mean KH are important for the purpose of upscaling (section 1.4), KG = exp (µ), (1.15)  1 KA = expµ + 21 2= KG exp2 2, (1.16)  11 KH = expµ - 2 2= KG exp-2 2. (1.17) The second moment of K(x) is given by   Var(K) = exp(2) - 1exp2µ + 12 2. The correlation structure of conductivity in space is captured by the covariance model. It can be expressed using Y ’s covariance model  Cov [K(x + s),K(x)]=exp 2µ + 2 + CVY (s). The statistical description of spatial correlation is based on the assumed covariance model. A bunch of models exist, which all exhibit properties of spatial correlation that can be observed in practice. Finding an appropriate model of correlation structure from field data includes an estimation procedure, see e.g Samper and Neuman [1989a, b]. In general, spatially distributed measurements are log-transformed and interpreted to gain an experimental semivariogram or a covariance function, respectively. The data is then fitted to the covariance model. At this point different models might lead to equivalent good results describing the correlation structure. The exponential, Gaussian and spherical models are typical representatives for covariance models used in practice [Rubin, 2003]. They differ in the near-origin behavior and the rate of decay of correlation as function of distance. The exponential model is indicative of sharp transition occurring between blocks of different values and the spatial correlation decays at a very low rate with large distances. The Gaussian covariance model describes a gradual transition between blocks of different conductivity values and shows a faster decay of correlation. The spherical model describes sharp transitions and discontinuities between zones of different conductivity 1. Principles of Fluid Flow and Salt Transport in Heterogeneous Porous Media values. In practice it is often very difficult to determine which of the models fit best to the observed data, e.g. Turcke and Kueper [1996] . Riva and Willmann [2009] showed that mean values of quantities like the hydraulic head are not considerably influenced by the choice of one of these three variogram types. We focus on a 3D-Gaussian model with the auto-covariance function CovG(s)= 2G(s), where G(s) is the correlation function X G(s) = exp -(si/`i)2 . (1.18) i The model relies on stationarity and hence can be expressed as function of the separation distance s =(s1,s2,s3). The Gaussian model refers to a finite correlation structure, which is expressed by the correlation length . =(`1,`2,`3) > 0. It characterizes the degree of continuity in every direction. It is most likely to assume that the correlation length is different for different directions. In case of equality `1 = `2 = `3 the medium is said to be isotropic. Otherwise it is called anisotropic, where the anisotropy between the two directions i and j is expressed by the ratio e = `i/`j. Another quantity commonly used to describe the spatial persistence of K(x) is the integral scale defined as Ii =1/2R8 C(si)dsi, in ith direction. The larger the integral scale, the 0 longer is the spatial persistence of the correlation. The integral scale is proportional to the correlation length and can be transformed according to the model of choice. For the Gaussian model the relation is given by Ii = v /2`i. These mentioned models all assume a finite correlation scale continuously distributed over space. However, geological media exhibit correlation on multiple scales depending on the problem of consideration and observation [Gelhar, 1986; Dagan, 1986]. A fractal model allows to account for an evolving heterogeneity structure. Neuman et al. [2008] showed that, even though sample data appear to fit either Gaussian or Exponential variograms, they could be represented equally well in terms of a truncated power law variogram model. The semivariogram of the power law model is given by SVPL(r)= C0r 2H , (1.19) with Hurst coefficient 0 0.5 the probability of observing increments to the same sign over large distances is high and fields becoming smooth in appearance. For H< 0.5 larger small-scale variability can be observed and the appearance is more rugged. The unbounded growth of variability is a theoretical hypothesis which cannot be captured in numerical studies. For generating fields of hydraulic conductivity numerically a truncated power law is used, where a cutoff length is chosen larger than the domain size. 1.4. Upscaling 1.3.3. Quantifying Variables of Flow and Transport in Stochastic Media Two approaches are commonly used to solve a problem in stochastic subsurface hydrology: the analytical and the numerical approach. The idea of the analytical approach is to derive explicit expressions for effective quantities from stochastic differential equations. The numerical approach uses random field generators to produce hydraulic conductivity fields. They serve as input to numerical flow and transport models which provide numerical solutions for the flow field and/or concentration patterns in heterogeneous media. Since hydraulic parameters may vary over distance and time, variables of flow and transport like hydraulic head or the concentration distribution of a tracer may be spatially dependent, too. Introducing K(x) as SRF into Darcy law, Eq.(1.2), and the mass conservation equation, in Eq.(1.5), results in stochastic PDE’s. Hence, flow variables are also SRF. It is often aimed to express these SRF with the statistical parameters of K(x), which are KG, 2 and f for a Gaussian covariance model. In general, the hydraulic head is less variable than the log conductivity due to physical constrains, given by the PDE’s [Rubin, 2003]. Using the Monte Carlo approach to solve the flow equation with heterogeneous conductivity distributions numerically is conceptually simple and needs no particular assumption. It allows the analysis of complex scenarios which are not accessible to analytical methods. The main steps of such a procedure are (i) assuming a probability distribution of the known input variable, which is in general the hydraulic conductivity K(x); (ii) generating realizations synthetically; (iii) solving the flow equation with appropriate boundary conditions for every realization; (iv) computing the statistics of the model output variables. An advantage of this approach is the flexibility in the problem description. A disadvantage is that the calculation of huge ensembles of realizations might be necessary, which can be quite consuming in time and computational power. Furthermore, it is often difficult to predict how many realizations and thus simulations are necessary to achieve convergence of solutions. 1.4. Upscaling One of the questions arising in groundwater flow problems is how to treat heterogeneity in aquifer models. In order to reduce complexity it is important to find representative parameters, which are capable of reproducing some average behavior of the system by upscaling. Upscaling basically aims to connect scales, either to connect conceptual model scale to observation scale or to connect different observation scales. The emergence of upscaling in hydrogeology is directly linked to the application of geostatistical methods on subsurface flow and transport problems. Matheron [1967] was the first who described heterogeneously distributed hydraulic conductivity as SRF and directly asked for representative values and bounds. Since then a huge amount of literature emerged, concerning upscaling of hydraulic conductivity. For a detailed review see e.g. Dagan [1989]; Renard and de Marsily [1997]; Sánchez-Vila et al. [2006]. In the context of stochastic subsurface hydrology, numerical models can serve for validating 1. Principles of Fluid Flow and Salt Transport in Heterogeneous Porous Media theoretical upscaling approaches. It is often very complicated to define upscaled hydraulic parameters for complex flow regimes. Due to the lack of data it is very difficult to validate averaged values obtained from upscaling approaches directly. Numerical models can be used to fill the gap of missing measurements. Based on the physical process equations we can perform numerical experiments and generate virtual data. Results from upscaling theory can be validated by comparing analytical solutions with numerical experiment data. The first following this strategy was Freeze [1975]. We discuss the combination of theoretical framework and numerical methods to find upscaled descriptions of hydraulic conductivity in the flow conditions of uniform flow and well flow in detail in the sections 2.3 and 4, respectively. 1.4.1. Scale Hierarchy Heterogeneity in hydraulic properties of porous media can be observed on multiple scales. The characterization of the spatial correlation scale of the property of interest coincides with the definition of spatial scales of flow domains. Dagan [1986] defined three major scales: the laboratory, the local, and the regional scale with corresponding correlation on pore scale, correlation of hydraulic conductivity, and correlation of transmissivity. Heterogeneity on pore scale is not considered here. For flow on local scale the hydraulic conductivity K(x)= K(x, y, z) is the variable of interest. The vertical component of the flow has to be taken into account, since the flow patterns are strongly influenced by the three dimensional nature of the porous media. Furthermore, anisotropy effects have to be considered. For flow on regional scale the relation between the horizontal extent and the aquifer thickness becomes very small. Thus, it is assumed that vertical flow can be neglected and the flow field can be seen as horizontal. The porous media can be described by the two dimensional transmissivity T (x), which is the vertical average of the hydraulic conductivity K(x)= K(x, y, z) over the domain thickness L, ZL T (x, y)=K(x, y, z)dz. 0 1.4.2. Effective and Equivalent Conductivity Two main approaches have evolved in literature to find representative parameters. The first is the approach of effective quantities which are gained by performing ensemble averages. The second approach uses spatial averaging to find equivalent parameters. However, determining upscaled parameters for a certain problem of consideration can be quite challenging. Effective values might not exist or spatial averages can strongly depend on the scale of the averaging volume. Effective hydraulic conductivity is defined as the value Kef (or Keff) which satisfies Darcy’s law hq(x)) = -Kefhrh(x)iwith ensemble averaged quantities h.iof flux q and hydraulic head h. The effective conductivity relates expected values of specific discharge and head gradient. Kef is defined as a function of aquifer’s material properties and must be valid throughout the entire 1.4. Upscaling domain; it is a characteristic property of the medium but must not be influenced by flow conditions [Rubin, 2003]. Closed-form results for an effective conductivity Kef are available in steady uniform flow for cases of low variability and simple flow regimes, discussed in more detail in section 1.4.3. In complex flow configurations or for large variances the effective conductivity can be bracketed by bounds of arithmetic and harmonic mean KA and KH. These are given in Eq.(1.16) and Eq.(1.17) for log-normally distributed hydraulic conductivity fields [Matheron, 1967; Sánchez- Vila et al., 2006]. Non-uniform flow regimes like in pumping tests lead to Kef depending on space coordinates and/or boundary conditions. Thus, the common definition of effective conductivity is not appropriate to describe representative conductivity in these flow regimes. Another approach to define an upscaled conductivity refers to volume averaging. The equivalent ¯¯r-= Kq eq h with the spatial average of the Equivalent parameters are associated with conductivity Keq is defined to fulfill Darcy’s law R and of the head gradient h¯ . flux q = a particular geometry and boundary conditions and can thus be defined for non-uniform flow conditions. In case of ergodicity, ensemble averages can be replaced by spatial averages. Hence, effective and equivalent conductivity are identical. In stationary conductivity fields and when flow is induced by a linear pressure gradient, ergodicity can be assumed throughout the entire domain. The mean flux is equal to the spatial average of the flux [Rubin, 2003] and thus Kef = Keq. For non-uniform flow this relation is not valid, due to the ergodicity breakdown. Beyond the two approaches of effective and equivalent parameters a lot of work has been devoted to finding representative descriptions of conductivity for more complex flow regimes, as discussed in detail e.g. in Sánchez-Vila et al. [2006]. ¯ 1.4.3. Upscaled Conductivity for Uniform Flow Steady state parallel flow in stationary media is one of the particular situations where effective conductivity, denoted by Kefu, exists. The theoretical approach provides analytical solutions for the effective conductivity of a log-normally distributed isotropic conductivity K(x) in dependence on the dimensionality d of the considered flow field, 11 Kefu = KG exp 22 - , (1.20) d with the statistical quantities geometric mean KG (Eq.(1.15)) and variance 2 of K(x). For the more complex solution in anisotropic media, see section 4.2.2 or [Gelhar and Axness, 1983; Dagan, 1989]. For all dimensions Kefu is bounded between KA and KH, regardless of the assumed correlation (1d) structure. In one dimension the effective conductivity is the smallest with K= KH = efu  V q dV 1 V KG exp- 2, due to the fact that in one direction a serial flow pattern dominates, where the harmonic mean is the average of choice. In two and even stronger in three dimensions, flow 1. Principles of Fluid Flow and Salt Transport in Heterogeneous Porous Media (2d) can circumvent areas of low conductivity, thus the effective values are given by K= KG efu  (3d) = KG exp 16 2 . (or transmissivity Tefu = TG, respectively) and K efu In section 2.3 we show that the analytical description of Kefu in Eq. (1.20) can be reproduced by the effective hydraulic conductivity calculated from an ensemble of simulations with uniform flow through two-dimensional heterogeneous log-normally distributed conductivity fields. 1.4.4. Upscaled Conductivity for Well Flow Most field methods, which rely on well flow scenarios to obtain hydraulic conductivity values, assume homogeneity [Kruseman and de Ridder, 2000]. In homogeneous media well flow can be described by Thiem’s solution [Thiem, 1906] for steady state h(r)= - Qw ln r + hR (1.21) 2LK R or Theis’ solution [Theis, 1935] for transient flow, where h(r) is the hydraulic head, depending on the radial distance r from the well. Qw is the discharge at the well, L is the aquifer thickness, K is the homogeneous conductivity, and hR = h(R) is the known hydraulic head at a certain distance R from the well. The analytical solution of Thiem (1.21) can be used to benchmark numerical codes for convergent flow regimes, see section 2.2. Well flow in heterogeneous media reveals a much more complex behavior than in homogeneous media. The constitutive equation, which is the average of Darcy’s law, has a non-local structure for non-uniform flow and it is in general not possible to find a single upscaled conductivity [Matheron, 1967; Indelman and Abramovich, 1994; Indelman et al., 1996]. The steep pressure gradients which occur at the well, lead to a breakdown of the ergodicity assumption. An effective conductivity value Kef, which is valid throughout the entire range, does not exist and thus the definition of effective parameters is not applicable to convergent flow. Although a single representative conductivity for well flow does not exist, the asymptotic behavior is well known [Indelman and Abramovich, 1994]. Flow can be assumed to be uniform  in average far from the well, thus the representative value is Kefu = KG exp 16 2 in three dimensions and Tefu = TG in two dimensions. The near-well representative value Kwell or Twell is stronger impacted by heterogeneity due to the high pressure gradients. The value depends on the boundary condition assumed at the well. It is either the arithmetic mean KA or TA for a constant head or the harmonic mean KH or TH for a constant flux boundary condition, e.g. [Shvidler, 1966; Dagan, 1989; Indelman et al., 1996; Indelman and Dagan, 2004]. A detailed discussion on that is given in section 4.2.2. Since ensemble averaging is not applicable, the concept of spatial averaging can be used to obtain an equivalent conductivity, which depends on the distance to the well. Indelman et al. [1996] defined the equivalent conductivity Keq(r) based on Thiem’s formula (1.21) as a spatial 1.4. Upscaling average over a circle of radius r around the well in the sense of Matheron [1967], Qw(ln r - ln rw) Keq(r)= 2L(h(r) - h(rw)), where Qw is the flow rate, rw is the radius of the well, h(rw) is the head at the well and L is the aquifer thickness. Keq(r) can be interpreted as the conductivity of a cylinder of radius r in homogeneous media which conveys the same discharge as the actual medium for the same mean head difference. In contrast to Kefu in uniform flow Keq(r) depends on the flow geometry. Firmani et al. [2006] showed by the use of numerical simulations that Keq(r) is appropriate to describe the conductivity of well flow only for small variances and high anisotropy ratios. A large amount of work is focused on finding upscaled descriptions of the conductivity and transmissivity in well flow regimes expressed by a radial dependency. Some studies are based on the definition of Keq(r) or Teq(r), respectively; e.g. Fiori et al. [1998]; Dagan and Lessoff [2007]. But also other approaches have emerged, e.g. Neuman and Orr [1993]; Sánchez-Vila [1997]; Sánchez-Vila et al. [1999]; Guadagnini et al. [2003]; Neuman et al. [2004]. 1.4.5. Well Flow Adapted Filtering Schneider and Attinger [2008] presented a promising alternative approach of upscaling well flow in two dimensions. They presented a radially dependent upscaled transmissivity T CG(r), which allows a transition from near field transmissivity TH to far field value Tefu = TG controlled by the correlation length `, . 1 2 . T CG (r)= TG exp - , (1.22) H 2 (1 + 2r2/`2) where r is the radial distance from the well, TG is the geometric mean, 2 is the variance and . is the correlation length of the log-normally distributed transmissivity T (x). . is a factor of proportionality. The formula of T CG(r) in Eq(1.22) was derived making use of the upscaling procedure Coarse H Graining, introduced for uniform flow by Attinger [2003]. The procedure of Coarse Graining and its application on two-dimensional well flow can be found in Schneider and Attinger [2008]. We discuss it for three-dimensional well flow in section 4.3. The fundamental idea is to perform a spatial filtering on the flow equation which is appropriate to the non-uniform character of a pumping test. The filter is proportional to the radial distance from the well. Hence near the well the filter length is very small, so nearly no filtering is applied and the heterogeneity of the local transmissivities is still resolved. Far away from the well the filter volumes are very large and the local heterogeneous transmissivity values are replaced by the effective value Tefu. The analytical expression of T CG(r) in Eq.(1.22) allows a direct calculation of the corresponding H upscaled hydraulic head hewf(2d)(r) by solving the groundwater flow equation, as done in the Appendix A.4. The solution hewf(2d)(r), called the effective well flow head, provides an effective 1. Principles of Fluid Flow and Salt Transport in Heterogeneous Porous Media description of the hydraulic head for large-scale pumping tests in a two-dimensional flow regime  hewf(2d)(r)= - Qw ..(z(r)) - ..(z(R)) - e  22 ..z(r) - 2 +e  22 ..z(R) - 2 + hR , 4TG22 (1.23) where r is the radial distance from the well, Qw is the pumping rate, TG is the geometric mean, and 2 is the variance of the log-transmissivity. ..(z) is the Theis function or exponential integral R z exp(z0) 01 function ..(z)=0dz. z(r) is an abbreviation for z(r)= 2 (1+(r/`)2), likewise z(R)= -8 z2 2 1 (1+(R/`)2), with f being the correlation length and . being the factor of proportionality. R is an arbitrary distance from the well, where the hydraulic head h(R)= hR is known. Numerical simulations of pumping tests in a two-dimensional flow regime showed that ensemble averaged drawdowns can be reproduced very well by making use of T CG(r) as shown by H Schneider and Attinger [2008]. 2. Numerical Benchmarks 2.1. Introduction Complicated hydrogeological problems can be analyzed using modern numerical codes. Major challenges for scientists are the understanding of the coupled processes, the implementation of algorithms, and the integration of the available experimental data. To cope with the problems, interdisciplinary cooperation and an interactive development procedure are required. Benchmarking is a common procedure for model verification. Complex numerical models rely on a properly and efficiently implemented code. Therefore, the numerical software is tested on examples. In case analytical solutions exist, benchmarking comprises of the comparison with the numerical solution. Another approach is to compare simulated results with observational data. In these cases caution has to be given to the fact that measurements can be biased. For models without analytical solutions or observational evidence a cross comparison with other numerical codes is often done. For the simulation software OpenGeoSys an extensive benchmarking documentation is given in Kolditz et al. [2012a]. The scope of this section is to introduce four specific benchmarks that refer to the complex problems discussed in parts II and III. On the one hand, it serves as code testing. On the other hand, we outline important features of numerical simulations dealing with heterogeneous parameter distributions and density dependent flow. The four benchmarks can be understood as examples for the development of simple hydrogeological models with reference to the steps of model evolution, as described in section 1.1. We introduce the conceptual model and state the process equations of relevance which are discussed in detail in section 1.2. We further adapt the problems to the software OpenGeoSys [Kolditz et al., 2012b]. This includes defining the model domain with a spatial grid resolution, fixing a temporal resolution, defining initial and boundary conditions as well as specifying physical parameters. In addition, the discussion of each of the four benchmarks includes a short literature review. 2.2. Pumping Test in Homogeneous Porous Media Pumping tests are a widely used tool to identify hydraulic properties. In a pumping test water is extracted from a single location which causes a radial flow character. From the mathematical point of view this point source creates a singularity at the pumping well, which complicates the investigation of the flow equation in particular in heterogeneous porous media. For homogeneous media an analytical solution exists, introduced by Thiem [1906]. We benchmark the numerical solver OpenGeoSys for convergent flow by comparing a simulated pumping test drawdown with Thiem’s solution. 2. Numerical Benchmarks Figure 1: Radial refinement of the numerical mesh in the area of 1m around the pumping well. 2.2.1. Thiem’s Solution We give a short derivation of Thiem’s solution. The starting point is the equation for saturated groundwater flow, which is a combination of the continuity equation of flow (1.1) and Darcy’s law (1.3), @h(x,t) S -r· (K(x)rh(x,t)) = Q(x,t) , (2.1) @t where S is the storage capacity, h(x,t) is the hydraulic head in space x =(x1,x2) and time t, K(x) is the hydraulic conductivity and Q(x,t) represents a source or sink term. For a steady state pumping test the time derivative is zero. Furthermore, we assume homogeneity in hydraulic conductivity, thus K(x)= K. Transforming Eq.(2.1) from Cartesian to polar coordinates (x1,x2) . (r, ), with radius r and angular coordinate , results in 1dh(r)d2h(r) 0= K + K 2 . (2.2) r dr dr According to the nature of radial flow the solution h(r) does not depend on the angular coordinate . In Eq.(2.2) the sink term is not included, which is instead used as boundary condition. A constant discharge Qw at the well corresponds to a Neumann boundary condition and is calculated as the flux over the surface w of the pumping well, I Qw =q(rw)d w = -2rwLKh0(rw), w where rw is the radius of the well, L is the aquifer thickness, and q(r) is the flux, which is given by Darcy’s law in radial coordinates q(r)= -Kh0(r)er with er being the unit vector in radial direction. The second boundary condition for the ODE (2.2) is given by a fixed hydraulic head h(R) at the radial distance R. 2.3. Uniform Flow in Heterogeneous Porous Media Figure 2: Drawdowns of simulated pumping test (dots) and analytical Thiem’s solution (solid line) in dependence on the radial distance r from the well. The solution of the ODE (2.2) is given by h(r)= C1 ln r + C2. The integration constants C1 and C2 are determined using the previously discussed boundary conditions, resulting in C1 = dh(r) Qw Qw r = - and C2 = ln R + h(R). Thus Thiem’s solution reads dr 2LK 2LK h(r)= - Qw ln r + h(R) . (2.3) 2LK R 2.2.2. Simulation Results We carried out simulations on a two dimensional numerical mesh ranging from -6.4m to 6.4m in both directions. We established a radial refinement at the well in the range of -1m to 1m as visualized in Figure 1. It provides high resolution of the hydraulic head drawdown around the well. The well is not included as a point source but as a hollow cylinder with radius rw =0.01m. We chose a constant head of zero at the outer radial distance R =6.4m as boundary conditions. At the well we used a constant total pumping rate of Qw = -1e-3m3/s. Figure 2 shows the comparison of the drawdowns from a simulated pumping test and the analytical Thiem’s solution, Eq. (2.3). Theoretical and numerical drawdowns are in very good agreement. Especially at the well the refined mesh provides good resolution. Thus, the simulation software OpenGeoSys can reproduce very well the radial dependent groundwater flow with steep gradients at a pumping well. 2.3. Uniform Flow in Heterogeneous Porous Media The existence of an analytical expression for the effective hydraulic conductivity Kefu makes the uniform flow regime favorable for numerical benchmarking. We performed simulations in heterogeneous porous media with a uniform flow field in a simple domain design. It allows the 2. Numerical Benchmarks application of the Monte Carlo approach to determine the statistics of flow in stochastic media, as described in section 1.3.3. From the simulated hydraulic head distribution a simulated effective conductivity Kˆ efu can be inferred using Darcy’s law and averaging strategies. The benchmark focuses on the comparison of the analytical expression of Kefu with the simulated effective conductivity Kˆ efu. Analytical and Simulated Effective Hydraulic Conductivity The evaluation of an effective conductivity for uniform flow in heterogeneous porous media has been subject of several studies since the beginning of geostatistics in hydrogeology. Matheron [1967] started by obtaining the upper and lower bounds, later Gelhar and Axness [1983]; Dagan [1989] presented a general analytical expression, called Kefu. The expression for all dimensions d =1, 2, 3 can be written according to Eq. (1.20) as 21 1 Kefu = KG exp 2 - , d with KG being the geometric mean and 2 the variance of the log-normal distributed conductivity K(x). Several authors addressed the issue of comparing the analytical expression for Kefu with effective mean values obtained from simulations. Although the analytical expression of Kefu was derived for unbounded domains, it can be assumed that it is valid in bounded domains for a sufficiently large domain extent, discussed in detail in Sánchez-Vila et al. [2006]. Freeze [1975] was the first who performed numerical investigations of one dimensional flow in heterogeneous media. Follin [1992] gave a numerical confirmation of Kˆ efu = Kefu in two dimensions (up to 2 = 16) and Dykaar and Kitanidis [1992b] in three dimensions (up to 2 = 6). Dykaar and Kitanidis [1992b] also investigated the relation of domain extent and convergence of the effective conductivity toward an asymptotic value. Their results confirm that an analysis with Kefu is valid in case of large domain extent with respect to the corresponding directional correlation length. To obtain an effective conductivity Kˆ efu from simulations, we averaged the simulated hydraulic Rˆ head hˆ in the direction perpendicular to the flow direction, hˆ (x)=h(x, y)dy. In the second step, we inserted hˆ (x) into an inverted Darcy’s law, Eq.(1.3), Z x2 (Lx - x) Kˆ efu =Qw dx, x1 hˆ (x) where Qw denotes the flow rate and Lx the horizontal domain extent. The integration limits x1 and x2 were chosen sufficiently far from the boundaries. Numerical Simulation Setup We carried out simulations in a two-dimensional domain with a spatial extension of Lx = 200 m and Ly = 50 m and a regular grid spacing of 1m. We established a constant flow rate of Q = 2.3. Uniform Flow in Heterogeneous Porous Media Figure 3: Realization of a generated normally distributed log-conductivity field Y = ln K(x) with Gaussian correlation function. Light patches indicate high conductivity, whereas dark patches indicate low conductivity. Statistical parameters: KG = 1e-4m/s, 2 = 1, f = 5 m, and e = 1. 5·e-7m3/s at the left boundary (x = 0m) as inflow. The outflow area is at the right side (x = 200 m), with a constant head of h(x = 200) = 0 m. The top and bottom boundaries are impervious to the flow. The setting resembles a section of a confined aquifer with a mean head gradient of 1m over a distance of 200m. We generated an ensemble of 20 heterogeneously distributed conductivity fields. Fields are of log-normal distribution with a Gaussian correlation structure, as defined in Eq. (1.14) and Eq.(1.18). The statistical parameter are the following: mean conductivity of KG = 1e-4m/s, variance of 2 = 1, correlation length of f = 5 m, and anisotropy ratio of e = 1. We visualized a snapshot of one realization in Figure 3. One correlation length is resolved over five elements of the numerical grid. The relation of domain size to correlation length is Lx/f = 20 which is sufficient for Gaussian correlated fields to neglect boundary effects according to Dykaar and Kitanidis [1992b]. Results We plotted the simulated drawdowns of 20 realizations and the ensemble mean drawdown in comparison to the drawdown in homogeneous media with Kefu = KG as substitute value in Figure 4. The effective conductivities Kˆ efu of single realizations deviate up to 10% from the analytical solution Kefu = KG. The ensemble mean of the 20 simulation runs deviates only less than 1% from KG. In general, it was possible to reproduce the effective conductivity for uniform flow Kefu with numerical simulations of flow in heterogeneous media quite well. However, the domain size is too small to guarantee convergence to Kefu within a single realization. Convergence to the analytical solution can be found by taking the mean of several simulations runs. These results underline the fact that the sample size and the number of realizations is crucial in simulations with heterogeneous conductivity distributions. 2. Numerical Benchmarks Figure 4: Simulated drawdowns of hydraulic head: drawdowns of 20 simulation runs (dotted lines), ensemble mean drawdown (dashed line) and drawdown in homogeneous media with Kefu = KG as homogeneous substitute value (solid line). 2.4. The Elder Benchmark The Elder benchmark is one of the commonly used benchmarks to verify density-dependent flow in hydrogeological simulation software [Voss and Souza, 1987; Oldenburg and Pruess, 1995; Kolditz et al., 1998; Frolkovic and Schepper, 2001; Diersch and Kolditz, 2002; Prasad and Simmons, 2003; Park and Aral, 2007]. The wide use of this benchmark is due to the fact that it is a problem of free convection where flow is purely driven by density differences. Background The benchmark is based on numerical studies of the "short heater problem" introduced by Elder [1967a, b] in combination with experimental measurements of thermally driven convection in porous media. Voss and Souza [1987] adapted and reformulated the problem in terms of density dependent flow induced by salt concentration differences for code testing purposes. Resulting salt concentration distribution patterns of Elder [1967b] and Voss and Souza [1987] matched quite well. In subsequent numerical studies different steady state solutions, varying in number and development of salt plumes, have been observed e.g. Oldenburg and Pruess [1995]; Frolkovic and Schepper [2001]; Diersch and Kolditz [2002]; Park and Aral [2007]. At first, these discrepancies have been mainly attributed to mesh resolution [Frolkovic and Schepper, 2001; Diersch and Kolditz, 2002]. By making use of bifurcation analysis Johannsen [2003] showed that (at least) three stable and a further eight unstable steady state solutions co-exist. Different domain discretizations, numerical schemes and initial conditions can lead to different solutions 2.4. The Elder Benchmark Figure 5: Geometry and boundary conditions of the Elder benchmark for the left half domain. in number of down-welling plumes, shown by Park and Aral [2007]. However, the existence of multiple steady state solutions is an intrinsic characteristic of the Elder problem. The lack of a unique steady state solution makes the Elder problem questionable for benchmarking purposes. By making use of a pseudo-spectral approach, van Reeuwijk et al. [2009] confirmed the results of Johannsen [2003] and underlined the fact that bifurcation depends on the Rayleigh number as defined in Eq. (1.11). According to their analysis only one steady state solution exists for lower Rayleigh numbers than originally used in the Elder benchmark. Consequently, they established a low Rayleigh number Elder benchmark with slightly changed setting. We adopting that idea and present results for both, the classical Elder benchmark as well as the low Rayleigh number Elder benchmark. Problem Definition The Elder domain is a vertical cross section of x = 600 m length and z = 150 m height (Figure 5). Half of the surface contains a salt source with a constant concentration of C =1. The bottom has a fixed concentration of zero. The hydraulic head is constant in the entire domain and all boundaries are impervious to flow and mass transport. The hydraulic parameters are listed in Table 2. Voss and Souza [1987] proposed to use these values to obtain a Rayleigh number of Ra = 400, following the lines of the "short heater problem" of Elder [1967b]. Following the lines of van Reeuwijk et al. [2009] the low Rayleigh number Elder benchmark remains with the geometry and the parameter setting of the classical Elder benchmark with the exceptions of: (i) a reduced value of maximal density of 1030 kg/m3 instead of 1200 kg/m3, resulting in a Rayleigh number of Ra = 60;(ii) a change in boundary condition at the surface: the concentration it fixed to zero at x =0 - 150 m and x = 450 - 600m; and (iii) a coarser time stepping. Parameters different for both simulation settings are listed in Table 3. 2. Numerical Benchmarks symbol [unit] quantity value f [-] porosity 0.1 . [m2] permeability 4.845e-13 C [-] relative concentration difference 1 0 [kgm-3] fresh water density 1000 . [kgm-1 s-1] viscosity 0.001 Dm [m2 s-1] molecular diffusion coefficient 3.565 e-6 l, t [m] solute dispersivities 0 Table 2: Parameters valid for classical and low Rayleigh number Elder benchmark. symbol [unit] quantity Ra = 400 Ra = 60 C [-] volumetric expansion 0.2 0.03 t [s] time stepping 1314000 7884000 vmax [m/s] maximal velocity 2e-6 3.0e-7 Peg grid Peclet number 1.3 0.2 Cr Courant number 1.1 1 Table 3: Parameters and numerical stability numbers for classical and the low Rayleigh number Elder benchmark. We reproduced the classical and the low Rayleigh number Elder benchmark with OpenGeoSys [Kolditz et al., 2012b]. We performed all simulations on the half domain, which is sufficient due to the symmetry of the problem. The mesh contains 128 · 64 identical square elements, corresponding to a mesh refinement level of 6, according to Frolkovic and Schepper [2001]; Diersch and Kolditz [2002]. The governing equations used to solve variable density flow consist of the three fundamental conservation equations: the continuity equation of flow, given in Eq. (1.1), the momentum equation (1.2), and the mass transport equation (1.5). The equations are coupled to the equation of fluid density (1.10) and the hydrodynamic dispersion (1.6). We used no numerical stabilization schemes. Calculating grid stability criteria like the grid Peclet number and Courant number, as e.g. defined in Kolditz [2012] result in values satisfying the criteria for numerical stability in spatial and temporal discretization (Table 3). Simulation Results We visualized the results of the concentration patterns in combination with the flow velocity field for both Elder Benchmarks in Figure 6. The first time step is at an early stage of plume evolution, the second step is at a later stage when simulations reached the steady state. For the Elder benchmark this is at t = 3 year and t = 30 years and for the low Rayleigh number Elder benchmark at t = 20 year and t = 200 years. The differences in evolution time corresponds to the difference in the Rayleigh number. 2.4. The Elder Benchmark Figure 6: Isolines of simulated salt concentration C = [0.1,...,0.9] of the classical and the low Rayleigh number Elder benchmark for two time steps: top left Ra=400, t=3years; down left Ra=400, t=30years; top right Ra=60, t=20years; down right Ra=60, t=200 years. The arrows mark the flow directions, their color the flow velocity. The evolving concentration pattern for the classical Elder benchmark shows several downwelling fingers at an early stage and finally one central down-welling finger at the steady state resembling the results of Frolkovic and Schepper [2001]; Diersch and Kolditz [2002]; Park and Aral [2007]. The results for the low Rayleigh number Elder benchmark correspond perfectly to the distribution of isolines presented in van Reeuwijk et al. [2009]. For a quantitative comparison we calculated the maximal penetration depth of several concentration isolines, listed in Table 4. For the penetration depth of the isoline of C =0.6 the value of the classical Elder benchmark corresponds perfectly to the calculated value given in Prasad and Simmons [2005]. The values for the low Rayleigh number Elder benchmark match the location of the isolines given in van Reeuwijk et al. [2009]. C PD (RaS = 400) PD (RaS = 60) 0.2 -145 m -138 m 0.4 -141 m -127 m 0.6 -136 m -113 m 0.8 -56.3 m -84.4 m Table 4: Maximal penetration depth PD[m] at steady state for different concentration isolines for the classical and the low Rayleigh number Elder benchmark. 2. Numerical Benchmarks Figure 7: Geometrical setting and boundary conditions of the Saltdome benchmark. 2.5. The Saltdome Problem The Saltdome benchmark reflects the idealized groundwater flow over a salt dome, taking significant variations of fluid density due to salt dissolution into account. It was originally formulated in the frame of the international HYDROCOIN groundwater modeling project (Case 5, Level 1 problem). The setting is in reference to the salt dome nuclear waste disposal site at Gorleben, Germany [Herbert et al., 1988; Oldenburg and Pruess, 1995; Konikow et al., 1997]. The model aims to identify how convective flow, which is induced by fluid density differences, interacts with the advective flow process initiated by pressure differences. We examine the Saltdome benchmark more closely since the idealized setting is of similar character as the cross section model of the Thuringian Basin, which we discuss in section 3. 2.5.1. Problem Definition The benchmark consists of three processes. Continuity equation of flow (1.1) and Darcy’s law (1.2) determine the groundwater flow. Eq. (1.5) describes the salt transport. The heat transport process requires the evolution of the energy balance equation in the porous media, Eq.(1.8). These processes are successfully implemented into the physical based simulation software OpenGeoSys [Kolditz et al., 2012b], which we used to evaluate the Saltdome benchmark. The model region is a two dimensional vertical cross section of 900m length and 300m height, representing a homogeneous and isotropic aquifer. The geometrical setting and the boundary conditions are visualized in Figure 7. We applied a linear pressure gradient of 1e5Pa at the top from left to right. We established a constant salt concentration boundary at the central third of the bottom, representing a salt source exposed to the aquifer. The boundary conditions of the heat transport process are given by a constant temperature of T =8 C at the surface and of T = 35 C at the bottom of the basin. It implies a constant temperature gradient 2.5. The Saltdome Problem Figure 8: Numerical mesh used for the Saltdome benchmark simulations. symbol [unit] quantity value f [-] porosity 0.2 . [mD] permeability 1000 . [kgm-1 s-1] viscosity 0.001 0 [kgm-3] fresh water density 1000 C [-] normalized salt concentration [0, 1] C [-] solute expansion coefficient 0.2 Dm [m2 s-1] molecular diffusion 5e-08 l, t [m] solute dispersivities 20;2 T [ C] temperature [8, 35] T [K-1] thermal expansion coefficient -3e-04 w [kgm2 s-3 K-1] thermal conductivity of water 0.6 s [kgm2 s-3 K-1] thermal conductivity of solid 2.0 cw [m2 s-2 K-1] thermal capacity of water 4200 cs [m2 s-2 K-1] thermal capacity of solid rock 1000 s [kgm-3] rock density 2700 Table 5: Parameters for the Saltdome benchmark in line with the problem definition of Herbert et al. [1988]; Kolditz et al. [1998]. of 30 C/km. The lateral and bottom boundaries are impervious to flow, salt, and heat. The initial concentration in the domain is C = 0. The initial temperature is T =8 C. The specific hydraulic parameters used in this application are given in Table 5. For numerical simulations, we use an unstructured grid consisting of 23581 triangular elements with gradual refinements in the area of the salt dome. The edge lengths are between 1m and 30 m. Figure 8 gives a visualization of the mesh. We adapt the time stepping to the choice of the mesh and the processes incorporated. Generally, a time step is in the range of one month. Simulations reached the steady state after 1000 years. We performed numerical simulations with two different approaches of transport. First, we assume conservative transport of salt and heat. We neglect their impact on fluid density by using a constant density model . = 0. The system is not coupled: mass and heat transport depend on fluid flow but not vice versa. Second, fluid density effects are taken into account. 2. Numerical Benchmarks Figure 9: Isolines of calculated salt concentration C = [0.1,..., 0.9] in the Saltdome benchmark domain after 1000 years simulation time for (a) the constant density model 0 and (b) for the concentration and temperature dependent density model (C, T ). The arrows mark the flow directions, the color of the arrows indicate the flow velocity. Differences in salt concentration and temperature affect the fluid density via a linear density model . = (C, T ), as described in Eq.(1.10). As a consequence, the fluid flow equation depends on the salt mass and heat transport resulting in a non-linear coupling of equations. 2.5.2. Results Simulation results of the salt concentration distribution for the constant fluid density model . = 0 are given in Figure 9a. A uniform regional groundwater flow regime develops from left to right due to the imposed pressure gradient. Mass and heat are transported conservatively with the groundwater according to the flow pattern. High concentrated solute (up to C =0.4) rises to the surface on the top right edge. We observe no salt at the left side of the salt dome area. The temperature pattern is similar, as it can be seen in Figure 10a. A stable temperature stratification is suppressed by inflowing cold water at the left side of the domain, whereas warmer water (up to T = 24 C) from the bottom region is transported to the surface in the discharge area at the right side. Figure 9b shows the simulated salt concentration distribution for the non-constant fluid density model . = (C, T ). In contrast to the results of the constant density model (Figure 9a), salt can be observed in at the left side of the salt dome. The entire bottom region of the domain is filled with highly concentrated brine. Though the concentration of brine which reaches the surface at the right edge is less than C =0.2. 2.5. The Saltdome Problem Figure 10: Calculated isolines of temperature T = [9, 12,..., 33] C in the Saltdome benchmark domain after 1 000 years simulation time for (a) the constant density model 0 and (b) for the concentration and temperature dependent density model (C, T ). The arrows mark the flow directions, the color of the arrows indicate the flow velocity. The regional groundwater flow is mainly present in the shallow basin region, whereas in deep basin regions the impact of the pressure gradient is significantly reduced. Flow occurs in vertical direction, visible at the flow field arrows. The temperature distribution resembles the salt distribution. The reduced impact of advective flow results in a smoother temperature stratification, as visible in Figure 10b. 2.5.3. Discussion A comparison of the simulation results for the different density models reveals that the evolving salt concentration and temperature patterns differ significantly from each other. The salt concentration is higher in the lower basin region for the fully coupled simulation with . = (C, T ), but less highly concentrated saline groundwater is transported to the surface in contrast to the uncoupled simulation with . = 0. The total amount of dissolved salt is higher for the concentration dependent density model, see Table 6. To emphasize on that the horizontal distribution of salt concentration at z = -200m depth C(x, z = -200m)forconstantandnon-constantdensitymodel isvisualizedinFigure11. It can be clearly seen that the slope of the salt plume is steeper for the constant density model, indicating a much stronger impact of advective flow in deep aquifer regions. In contrast to the advection dominated transport processes in the constant density simulations, 2. Numerical Benchmarks Figure 11: Salt concentration after 1000 years at z = -200m depth for the constant density model 0 (solid line) and the non-constant density model (C, T ) (crosses). 0 (C, T ) (C) (T ) C ¯0.097 0.229 0.224 0.086 Table 6: Mean concentration C ¯ for different density models, as defined in Eq.(1.12). the flow regime for salt and temperature dependent density is of mixed convection [Nield and Bejan, 1999]. The main groundwater flow is driven by a pressure gradient, but density effects can give rise to free convection at the bottom of the aquifer. This is supported by the high Rayleigh number of solute RaS ˜ 9600, calculated according to Eq.(1.11). The horizontal density differences provoke the formation of eddies, resulting in an increased velocity in vertical direction. Higher velocities in combination with dispersion lead to a stronger spread of salt and thus the effect of recirculation provokes stronger upconing. An increased density inhibits the upward transport of highly concentrated brine into shallow region by the groundwater flow. The evolving flow pattern represents a balance among buoyancy effects, advection and dispersion. Several studies [Herbert et al., 1988; Oldenburg and Pruess, 1995; Konikow et al., 1997; Kolditz et al., 1998; Diersch and Kolditz, 1998; Younes et al., 1999] investigated the effects of boundary and initial conditions, spatial discretization, dispersion and diffusion coefficients as well as other parameters for the Saltdome benchmark with concentration dependent fluid density. Of particular interest is the setting of diffusion and dispersion coefficients. Coarse numerical meshes require a certain numerical and/or physical diffusion or sufficiently large dispersion to prevent numerical errors and guarantees convergence. Large values of diffusion or dispersion results in stronger upconing of salt. The dependency of the amount of dissolved salt on the refinement of the mesh is due to the synergetic feedback effect of dispersive transport. It is more pronounced for coarser meshes since the numerically calculated advective flux depends on 2.5. The Saltdome Problem C [%] 0 (C, T ) (C) (C) 0 0 15.1 15.7 1.1 (C, T ) 15.1 0 0.6 15.2 (C) 15.7 0.6 0 15.7 (T ) 1.1 15.2 15.7 0 Table 7: Relative concentration difference C [%] between density models, as defined in Eq.(1.13). T [%] 0 (C, T ) (C) (T ) 0 0 18.8 19.5 1.9 (C, T ) 18.8 0 0.7 20.1 (C) 19.5 0.7 0 20.8 (T ) 1.9 20.1 20.8 0 Table 8: Relative temperature difference T [%] between density models, as defined in Eq.(1.13). the scale of discretization [Konikow et al., 1997]. Refined meshes allow for more realistic, lower values of diffusion and dispersion resulting in less salt spreading as shown by Younes et al. [1999]. In this line the discretization of the constant salt boundary condition at the bottom is of importance, since it regulates the amount of brine entering from the bottom source region. This in turn is responsible for the development of a gravity-driven clockwise circulation against the overall counterclockwise flow imposed by the pressure boundary conditions [Oldenburg and Pruess, 1995; Konikow et al., 1997; Kolditz et al., 1998]. The more brine is dissolved, the stronger is the recirculation effect. The impact of the temperature process on the evolving flow pattern is small. Simulation results for purely temperature dependent density . = (T ) (not visualized) and constant density . = 0 are very similar. The same is valid for simulation results for the variable density model without and with temperature . = (C) (not visualized) and . = (C, T ), respectively. These findings are underlined by the similar relative mean salt concentration and mean temperature, listed in Table 6 and the small relative difference in salt concentration and temperature distribution, listed in Tables 7 and 8. These findings are supported by the relatively low thermal Rayleigh number of RaT ˜ 170. It is above the critical value of Rac =4, which marks the onset of convection. However, advective flow suppresses the development of convection cells and thus the upward movement of fluid due to smaller density. This is in line with the findings of Diersch and Kolditz [1998]. Thus, temperature differences present in this setting are too small to significantly impact on the flow pattern. More recently, Holzbecher et al. [2010] investigated the sensitivity of double diffusive convection above a salt dome, in a setting comparable with the Saltdome benchmark. Holzbecher et al. [2010] confirmed that the formation of eddies is an effect of salinity differences. The smaller the buoyancy effect, the less pronounced is the upconing of salt. Moreover a stronger potential 2. Numerical Benchmarks flow reduces the tendency of eddy formation. Simulations accounting for temperature effects on fluid density underlined the fact that for gradients of 25 C/km, typically found in sedimentary basins in Germany, the effect of temperature on the evolving flow and salt pattern is marginal. Furthermore, Holzbecher et al. [2010] re-evaluated parameters aiming to use characteristic values most likely to be found in a real aquifer, confirming most values used in the Saltdome benchmark. Finally, we can state that the results for the Saltdome benchmark calculated with Open- GeoSys are in line with previous findings of Herbert et al. [1988]; Kolditz et al. [1998]; Diersch and Kolditz [1998]. Although the setting is highly idealized, the Saltdome benchmark is of general interest for studying strongly coupled groundwater flow and solute transport due to highly concentrated brine. The results show that flow patterns look significantly different when taking into account density effects. Even though the effects of mixed convection are different from a purely advective regime, the results point toward the fact the saline groundwater can be transported to shallow aquifer regions by an interplay of regional groundwater flow and dispersive flux in deep aquifer regions. Part II. Impact of Heterogeneity and Density Dependent Flow on Regional Scale 52 3. Mechanisms of Salt Transport in the Thuringian Basin 3.1. Introduction Saline groundwater, which comes close to or reaches the surface, is a phenomenon that can be observed in many places in the Thuringian Basin [Seidel, 2003]. However, it is not obvious why denser brine overlays lighter fresh water in this region. The hydrogeological processes, which cause the rising of saltwater plumes from deeper geological layers to the surface, are not yet fully understood. The aim of this numerical study is to investigate the mechanism of brine transport within the aquifers of shallow sedimentary basins in general, and of the Thuringian Basin in particular. Numerous mechanism can be responsible for fluid migration within sedimentary basins like gradients in hydraulic head driven by topography and fluid density variations, but also tectonic loading, seismogenic pumping and the production of diagenetic fluids [Person et al., 1996]. The relative importance and the degree of interaction of these driving forces on fluid flow varies, depending on the local hydraulic and tectonic conditions (e.g. permeability, porosity and geological structure). The role of specific driving mechanisms on fluid flow within different tectonic environments is examined in several studies, e.g. Magri et al. [2005, 2008] and Lampe and Person [2002]. The understanding of basin-scale flow and transport phenomena is important for geothermal applications, hydrocarbon migration and in particular for assessing long-term behavior of pollutants. Results from numerical simulations provide essential information about active processes in the subsurface on the basis of physical principles. Mathematical models complement field or laboratory-based investigations because they represent geologic processes that occur at very slow rates and over large length scales. Furthermore they take into account the simultaneity of flow, heat and mass transport and can consider coupling effects [Nield and Bejan, 1999; Simmons et al., 2001; Diersch and Kolditz, 2002]. Regional groundwater flow develops due to hydrostatic pressure gradients caused by differences in topography. Permeability is the property of porous media, which mainly determines the flow velocities and thus the groundwater flow pattern by Darcy’s law. However, the distribution of permeabilities is subject to a high degree of uncertainty. On the one hand, the amount of available measured data is very low. On the other hand, permeability is heterogeneously distributed in space. Variations over orders of magnitude are common [Aehnelt et al., 2011]. An additional factor of uncertainty are the permeabilities in fault zones. These regions often show very distinct hydraulic behavior compared to undisturbed basin compartments: sometimes they act as conduit, other times as barrier to groundwater flow [Lampe and Person, 2002; Bense and Person, 2006]. Due to the strong heterogeneity of permeability values, estimates of mean permeability on basin scale can hardly be representative for the entire area of interest. Generally, averages can only be given within ranges. Thus, large scale groundwater flow simulations using one set of homogeneous permeabilities for all sedimentary layers might not be appropriate to describe the fluid 3. Mechanisms of Salt Transport in the Thuringian Basin Figure 12: Geological surface map of Thuringia with focus on the Thuringian Basin (modified according to Thüringer Landesamt für Umwelt und Geologie (TLUG) [2002]). The bold lines indicate locations of cross sections, where A-A’ marks the location of the southwest-northeast cross section used for numerical modeling. dynamics satisfactorily. Furthermore, averaged permeability values for large sedimentary units do not take into account the effects of local heterogeneity on the flow pattern. In combination with density dependent flow, caused by salt concentration and/or temperature gradients, they can play an import role. Several studies like Simmons et al. [2001] point towards the fact that heterogeneity and disturbances in the geological layering can have quite diverse impacts on the onset and development of convection cells. Local differences in the permeability can cause the onset of convection. Contrariwise, heterogeneity can enhance mixing and thus lead to the breakdown of convection. Classical measures like the Rayleigh number are probably not appropriate to predict the onset of convection in heterogeneous porous media. The importance of temperature as source for density instabilities that can create free convection in deep permeable sedimentary rocks is well known, e.g. Elder [1967a]; Diersch and Kolditz [1998]; Nield and Bejan [1999]. Magri et al. [2005, 2008] showed for the North East Ger- man Basin that salinization is probably caused by the interaction of hydrostatic and thermal forces. The simulation results point out that the interaction between temperature induced deep thermohaline convection and topography driven flow in shallow aquifers can transport saline groundwater over large distances into near-surface regions. More recent studies by Kaiser et al. [2011] suggest that convection, which is triggered by temperature dependent fluid density variations, affects the flow pattern of the North East German Basin only locally. 3.2. Hydrogeological Characterization of the Study Area From Rayleigh theory it is known that the development of convection requires highly permeable sedimentary units with large vertical expansion and large differences in temperature [Nield and Bejan, 1999]. Comparing the flat and shallow geological setting of the Thuringian basin with the much deeper and strongly structured North East German Basin, it cannot be expected that thermal convection is the major driving force within the Thuringian basin. It is rather assumed that regional groundwater flow interacting with the fault structure is the major force driving deep saline water to shallow aquifers. To give numerical proof of the hypothesized fluid dynamics of the Thuringian basin, we performed transient coupled simulations of hydraulic, thermal and mass transport processes on a selected profile. The aim of the numerical study is to analyze the major mechanisms influencing the flow pattern: topography driven flow in a geological setting with faults, impact of hydraulic properties and heterogeneity, as well as fluid density differences caused by temperature and salt concentration gradients. 3.2. Hydrogeological Characterization of the Study Area 3.2.1. Geological Overview The geologically closed unit of the Thuringian Basin is located in the center of Germany. It stretches in an oval shape approximately 150km from northwest to southeast and 80km from northeast to southwest, see Figure 12. In the north, it borders on the Harz and Kyffhäuser mountains, in the west on the Eichsfeld Swell, and in the south on the Thuringian Forest mountains as well as Thuringian shale mountains [Hoppe, 1959; Jordan and Weder, 1995; Seidel, 2003; Kober, 2008]. Over much of its history the Thuringian Basin was part of the North German Basin. It became separated in the Late Cretaceous by the nearly 100km long Finne fault zone, which symbolizes the morphological and geological northeast border of the Thuringian Basin [Malz and Kley, 2012]. The recent tectonic structure of the Thuringian Basin developed mainly during the Late Cretaceous and the Cenozoic. The sedimentary succession of the Thuringian Basin started with the Upper Permian Zechstein. Polyphase deformations of the basements and of the sedimentary cover have created the complex tectonics of the Thuringian Basin from Triassic and latest Cretaceous time. Several northwest-southeast striking faults subdivide the synclinar structure [Seidel, 2003; Malz and Kley, 2012]. Their influence on the groundwater flow pattern can be significant on a local up to a regional scale. Depending on the size, location, direction and their petrophysical characteristics, they can act as barriers or preferential pathways. The basin is filled with Permian and Triassic sediments, while the overlying layers of Jurassic and Cretaceous age are mainly eroded. Tertiary and Quaternary deposits are restricted to recent river valley. The internal structure of the basin is strongly influenced by the Upper Permian Zechstein salt deposits. In the center of the basin the Zechstein is located in 1000 m depth, whereas the outcrops of the Zechstein deposits mark the edge of the basin [Seidel, 2003]. 3. Mechanisms of Salt Transport in the Thuringian Basin sedimentary unit sedimentary subunit abbreviation thickness hydrogeol. unit Trias Keuper middle km 0 - 220 m shallow aquifer lower ku 0 - 65 m Muschelkalk upper mo 45 - 75 m middle mm 45 - 115 m lower mu 85 - 120 m Buntsandstein upper so 100 - 190 m aquiclude middle sm 140 - 240 m deep aquiferlower Bernburg suB 150 - 170 m Calvörde suC 170 - 190 m Table 9: Diagrammatic plan of the sedimentary units used in the Thuringian Basin model with subunits, abbreviations, and thicknesses according to Gaupp et al. [1998] and Seidel [2003]. The rightmost column shows the hydrogeological classification which is based on the permeabilities of the sedimentary units, given in Table 10. 3.2.2. Hydrogeological Characterization of Sedimentary Units We reduced the hydrogeological characterization of the Thuringian Basin to the sedimentary units overlaying the Zechstein salt. It can be assumed that the thick sequence of Zechstein salt decouples mechanically and hydraulically the sub-salt Permian deposits from the Trias- sic sediments [Rödiger, 2005]. A diagrammatic plan with thickness and general hydrological characterization is given in Table 9. The dominating basin filling sediments are the Buntsandstein (Lower Triassic) and the Muschelkalk (Middle Triassic), which form the two main aquifer systems. Small areas of Keuper deposits operate as permeable shallow aquifers, mainly in the center of the Thuringian Basin. They have minor influence on the flow pattern due to their spatially limited occurrence. Younger remained deposits -from Tertiary and Quaternary -are not relevant for the large scale hydrogeology of the basin. To allow a differentiation of the lithology dependent hydraulic and thermal properties the Triassic units Buntsandstein and Muschelkalk were further divided into subunits (Table 9). Buntsandstein The Buntsandstein sediments in the Thuringian Basin have an overall thickness of 500- 720m [Puff et al., 2003]. They are structured into Upper, Middle and Lower Buntsandstein (so, sm and su). Due to significant differences in the hydraulic properties, the Lower Buntsandstein is further divided into Bernburg formation (suB) and Calvörde formation (suC), see Table 9. The Röt-formation (Upper Buntsandstein, so) contains large contents of clay, forming the regional most important aquiclude of the Thuringian Basin. The Röt separates the aquifer of Lower and Middle Buntsandstein from the Muschelkalk aquifer system [Jordan and Weder, 1995; Hecht, 2003; Rödiger, 2005], see Table 9. Only few values for measurements of permeabil 3.2. Hydrogeological Characterization of the Study Area 57 (1) [mD] (2) [mD] (3) [mD] (4) [mD] km 1 000 --- ku 1 000 --- mo 50 - 10 000 10 - 100 -- mm 50 - 10 000 10 - 10 000 -- mu 50 - 10 000 20 - 200 -2.5 so --10 0.025 - 2.5 sm 100 - 300 3 - 3 000 100 - 5 000 25 - 2 500 suB 100 - 300 3 - 3 000 100 - 5 000 25 - 250 suC 20 - 40 -100 - 5 000 2.5 Table 10: Permeability values . for the sedimentary units from literature: (1)Jordan and Weder [1995], (2)Hecht [2003] and Merz [1987], (3)Rödiger [2005], and (4)Aehnelt et al. [2011]. Transmissivities Tf from pumping tests in Jordan and Weder [1995]; Hecht [2003]; Merz [1987], and Rödiger [2005] are converted to permeability by . = Tf /(Lg). Permeability measurements by Aehnelt et al. [2011] from bore logs are multiplied by the factor 25 to account for effects of fractures according to Hauthal [1967]. ities can be found in literature, see Table 10. In his hydrogeological model of the Buntsandstein aquifer system in the Eastern Part of the Thuringian Basin, Rödiger [2005] used a conductivity of Kf (so)= 1e-7m/s. Bore-log permeability measurements [Aehnelt et al., 2011] point towards much smaller values of . =0.001 - 0.1mD, which corresponds to Kf = 1e-11-1e-9m/s. The facies of Lower and Middle Buntsandstein fluctuate between fluvial sandstones and lacustrine deposits [Gaupp et al., 1998]. The units of the Middle Buntsandstein (sm) and the upper part of the Lower Buntsandstein, the Bernburg formation (suB), can be assumed as one hydraulic unit [Jordan and Weder, 1995; Rödiger, 2005]. They form an aquifer with doubleporosity characteristic [Jordan and Weder, 1995]. The effects of the flow in the pore matrix is overlaid by the much faster reacting groundwater flow in the fissures and fractures. Comparisons of permeability measurements and pumping test results point toward a 20 - 30 times higher velocity in the fissures than in the pore matrix [Hauthal, 1967]. The fissured character of the sm/suB-aquifer is recognizable in large variations of permeability. Hecht [2003] stated a spread of permeabilities over four orders of magnitudes, depending on their position, from undisturbed rock formation up to strongly shattered sediments in crossings of fault zones. Recent investigations on Buntsandstein aquifer characteristics by Aehnelt et al. [2011] support these findings. In general, the range of permeability is very large, whereas values tend to decrease with depth. Typical values found in literature are listed in Table 10. The Calvörde formation (suC) of the Lower Buntsandstein has a mixed hydraulic character. The upper part, connecting to the Bernburg formation (suB) is acting as an aquiclude. Lower units, which connect to the Zechstein, consist of permeable deposits with aquifer characteristics [Jordan and Weder, 1995; Rödiger, 2005]. Thus, different mean values for the permeability emerge in literature, see Table 10. Permeability values of the suC are assumed to be smaller than those of the sm/suB-aquifer [Aehnelt et al., 2011; Jordan and Weder, 1995]. 3. Mechanisms of Salt Transport in the Thuringian Basin Muschelkalk The Muschelkalk is marked by a prevailing marine environment and is subdivided into three units which reach on average an overall thickness of about 250m [Gaupp et al., 1998]. The Lower Muschelkalk (mu) consists mainly of massive limestone with varying composition. The Middle Muschelkalk (mm) deposits are evaporites, such as dolomit marls, gypsum and even rock salt and dolomit limestone. The Upper Muschelkalk (mo) is characterized by limestones and marl [Gaupp et al., 1998]. The thicknesses of the subunits are listed in Table 9. The Muschelkalk forms a connected aquifer system, although it is marked by an alternating layering of aquifers and aquicludes [Hecht, 2003; Merz, 1987]. A closed thick aquiclude layer, like the Röt formation, is not present above or within the Muschelkalk. Since the Muschelkalk is mainly formed by limestone sediments, it is a Karst aquifer were flow takes place in faults and fissures. Depending on the degree of subrosion, large variations in permeability can be observed, which result in a strong heterogeneity of the aquifer [Merz, 1987; Hecht, 2003; Rödiger, 2005]. Higher permeabilities are often present where the rock formation shows a looser structure due to tectonic faults. Especially the low permeable limestone from the Muschelkalk can have increased permeable character in fault zones and change its character from aquiclude to aquifer [Jordan and Weder, 1995]. Mean permeabilities from literature are given in Table 10. 3.2.3. Hydrodynamic The hydrodynamic of the Thuringian Basin is dominated by the two large aquifer systems of the Muschelkalk and Buntsandstein. The aquiclude Röt separates both hydraulic units, allowing only little communication, except in fault zones or graben structures [Rödiger, 2005]. For the Buntsandstein areas the groundwater catchments coincide with the surface water catchments, giving that the groundwater is flowing mainly according to the morphological conditions [Jordan and Weder, 1995]. This can be assumed for the Muschelkalk aquifer system as well. The central part of the Thuringian Basin belongs to the Unstrut-catchment, whereas the eastern part belongs to the Saale-catchment. Areas of groundwater recharge are at the southwest, west and southeast margins of the basin. The main discharge takes place at the northeast margin, which coincides with the outflow direction of the rivers Unstrut and Saale. The study of Meincke [1967] points towards higher flow velocities at the recharge areas than in the basin center and in the discharge areas. Even a stagnating flow regime can be assumed for parts of the basin center. The different flow velocities indicate the existence of local spots of discharge. The general direction of groundwater flow is from southwest to northeast [Meincke, 1967]. The main fault zones striking from northwest to southeast can cause significant modifications of the flow pattern. A fault can act as barrier or preferential pathway depending on size, location, direction, and petrophysical characteristics. Thus, faults can connect or disconnect permeable layers and accordingly they can operate as areas of recharge or discharge [Meincke, 1967; Jordan and Weder, 1995; Hecht, 2003]. 3.3. Numerical Groundwater Flow Model Figure 13: Geological structure of the southwest-northeast cross section with fault zones and sedimentary units: km -Middle Keuper, ku -Lower Keuper, mo -Upper Muschelkalk, mm -Middle Muschelkalk, mu -Lower Muschelkalk, so -Upper Buntsandstein, sm -Middle Buntsandstein, suB -Lower Buntsandstein (Bernburg formation), suC -Lower Buntsandstein (Calvörde formation) . Superelevation factor of 10. Salinization of ground water due to subrosion of Zechstein salt was observed especially in the Saale-Valley [Hecht, 2003; Rödiger, 2005]. The effect of salinization is particularly strong in fault zones. 3.2.4. Temperature Distribution Hurtig and Meincke [1969] investigated the temperature distribution, thermal conductivities and heat flow in the Thuringian Basin. They showed that only small thermal anomalies are present in the Thuringian Basin. Unusual or disturbing heat sources are most likely not present. Observed heat flow anomalies are due to the relief and the rock composition of the basement. According to Meincke et al. [1967], inflowing meteoric fresh water influences the temperature distributionupto 400mdepth. Especiallyatthebasinmarginsacoolingeffectcanbeobserved. The thermal gradient amounts in average to 30 C/km. Higher gradients up to 33 C/km36 C/km were observed in elevated regions, like the Eichsfeld Swell. In flat areas like the basin center the gradients decrease to values of 21 C/km-27 C/km [Meincke et al., 1967]. The relatively low thermal gradients correspond to the fact that thermal water springs are not known in the Thuringian Basin [Hecht, 2003]. 3.3. Numerical Groundwater Flow Model 3.3.1. Cross Section Model We performed simulations on a southwest-northwest cross section of approximately 70km length and maximal 1000 m depth, shown in Figure 13. The profile ranges from the northern rim of the Thuringian Forest in the southeast, cutting through the region of Erfurt in the basin center, to the Finne fault in the northeast. The location of the cross section in the Thuringian Basin is depicted in Figure 12 as line A-A’. 3. Mechanisms of Salt Transport in the Thuringian Basin Figure 14: Detail of the numerical mesh in the environment of the Eichenberg-Gotha-Saalfeld fault zone. Sedimentary units are separated by bold black lines. Without superelevation. We chose the transect perpendicular to the northwest-southeast-striking faults. The choice is based on the one hand, on the general flow direction from southwest to northeast. On the other hand, we aim to investigate the effects faults have on the flow pattern, by incorporating the main fault structures, as given in Figure 13. The cross section model is based on the basin stratigraphy as described in Jordan and Weder [1995]; Gaupp et al. [1998]; Seidel [2003]; Rödiger [2005]. It contains the Triassic sediments of the Buntsandstein, Muschelkalk and Keuper as discussed in detail in section 3.2. The Zechstein deposits form a natural boundary in depth for the study area due to their impermeability to fluid flow. We did not incorporate the Zechstein layer into the cross section model. However, at the bottom of the basin, where the Triassic sediments are in contact with the Zechstein, dissolved salt can be transported by groundwater flow. We developed a work flow to gain a spatially discretized mesh, which is appropriate for numerical transport simulations from geostratigraphical input data. We generated a balanced crosssectionsmodelusingthesoftware2D- MoveonthebasisofgeologicalsurfacemapsGK 1 : 25000 of Thuringia with stratigraphical information from a digital underground model of Kober [2008]. We incorporated faults structures according to Malz and Kley [2012]. For generating the finite element mesh, we transformed and loaded the profile into the software Gocad. Subsequently we used Gocad and the interfaces described in Zehner [2011] to generate the numerical grid. A part of the numerical mesh in the environment of the Eichenberg-Gotha-Saalfeld fault zone is visualized in Figure 14. The unstructured finite element mesh consists of 121983 triangular elements. The average side length of a triangle is 35-40m. Themaximumsidelengthamounts to 50m. Locally refined triangles can have a side length of about 5m. The mesh satisfies the regularity conditions and the resolution allows to model variations in fluid density. We assume mesh convergence due to the fact that additional simulations, which we performed on finer meshes give the same results for all simulation scenarios. We incorporated discontinuities in the layering, which result from major fault structures, to the mesh. This results in a natural break up of the horizontal stratification into compartments of flat, undisturbed basin fill, interrupted by faults zones of regional importance (Figure 13). Partly the fault zones have a horizontal extension (Eichenberg-Gotha-Saalfeld fault zone, Erfurter fault zone and Finne fault zones), which allows for a separate assignment of parameters. 3.3. Numerical Groundwater Flow Model 3.3.2. Processes and Numerical Framework We carried out all simulations using the software OpenGeoSys. This finite element based simulator provides algorithms for solving coupled fluid flow, mass and heat transport processes in porous media. The software was successfully tested against a wide range of benchmarks for fluid, heat and mass transport by Kolditz et al. [2012a], and for density dependent flow in particular in section 2.4 and 2.5. Three processes have been established on the cross section model: a hydraulic process calculates the groundwater flow; a mass transport process gives the distribution of salt; and a heat transport process calculates the temperature distribution. Several partial differential equations describe these processes of flow and transport in saturated porous media: Darcy’s law given in Eq.(1.2), the continuity equation of flow (1.1), the equation of solute mass conservation (1.5) as well as the energy balance equation of the fluid and the porous media (1.8). The resulting system is fully implemented in OpenGeoSys. We conducted several simulation scenarios to investigate the impact of two main transport mechanisms. First, we performed simulations with different mean hydraulic permeability distributions. The parametrization of these scenarios is given in section 3.3.4. Thereby we used a constant density model of conservative mass and heat transport. Second, we investigated the role of density dependent flow. This required a coupling of all equations through a mass and temperature dependent fluid density model. In OpenGeoSys a linear relation for mass and temperature dependent fluid density is implemented, according to Eq.(1.10). 3.3.3. Boundary and Initial Conditions We set the fluid pressure at the surface boundary to a constant value of zero to establish a regional groundwater flow. This corresponds to a hydraulic head following the local topographical elevation. It induces a steady regional flow with recharge areas in the highlands and discharge areas in the lowlands. We implemented a mass transport in all simulations by fixing a constant salt concentration at the bottom of the basin model, where the Buntsandstein is in contact with the Zechstein salt. We fixed a zero mass concentration at the surface, representing the recharge with meteoric fresh water. Although this condition does not allow saline water to enter the surface, it is sufficient to localize zones where deep groundwater rises to near-surface levels. We applied a constant surface temperature of 8 C for the heat transport process. It corresponds to the current average annual temperature for the area under consideration. We fixed the temperature values at the lower boundaries corresponding to a linear vertical gradient of 30 C/km [Meincke et al., 1967]. Thus, the temperature at the bottom of the basin varies between 30 C and 38.4 C, depending on the depth and distance to the surface. We defined the lateral boundaries impermeable for fluid, heat and mass transport. The initial conditions are given by a zero pressure and pure fresh water distribution in the whole simulation 3. Mechanisms of Salt Transport in the Thuringian Basin symbol [unit] quantity value C - normalized salt concentration [0, 1] T [C] temperature [8, 38.4] T K-1 thermal expansion coefficient -0.0003 C - solute expansion coefficient 0.2 Dm [m2 s-1] molecular diffusion coefficient 1e-8 l, t [m] longitudinal and transversal dispersion 20; 2 0 [kgm-3] density of fresh water 1 000 . [kg m-1 s-1] viscosity 0.001 w [kgm2 s-3 K-1] thermal conductivity of water 0.6 cw [m2 s-2 K-1] thermal capacity of water 4 200 cs [m2 s-2 K-1] thermal capacity of solid rock 1 000 Table 11: Simulation parameters for the entire model domain. domain. The initial temperature was 8C. We performed all simulations with a computing time of 200000 years to assure that the system reaches steady state. 3.3.4. Parametrization Table 11 contains all model parameters for fluid, mass and heat transport processes, which are valid for the entire cross section model. Physical properties, which can be distinguished for the stratigraphical units, are listed in Table 12. We grouped the nine sedimentary layers to three main hydraulic units: a shallow aquifer consisting of Keuper and Muschelkalk (km/ ku/mo/mm/mu), the aquiclude unit Röt (so), and the deep aquifer formed by sediments from Middle and Lower Buntsandstein (sm/suB/suC), see also Table 9. Reasonable ranges and representative mean values of permeability are given in Table 12 for every sedimentary layer, based on the findings from literature as discussed in section 3.2 (Table 10). We considered each layer homogeneous and isotropic in density, thermal properties and porosity. We applied three reasonable mean permeabilities for every sedimentary layer: a lower (l), a reference (r) and a higher (h) permeability, given in Table 12. Having three permeability distributions for every of the three hydraulic units (shallow aquifer, aquiclude and deep aquifer) results in 27 possible simulation scenarios with different mean permeabilities. Additionally, we performed simulations with varied permeability values within the fault zones (EichenbergGotha- Saalfeld fault zone, Erfurter fault zone and Finne fault zone). While we kept the permeabilities in the undisturbed basin units constant at the reference value (r), we varied the permeability values within the faults over several orders of magnitude by multiplying with the factors 10, 100, 0.1 and 0.01. This is in line with findings of permeability values in fault zone in the Thuringian basin [Hecht, 2003; Aehnelt et al., 2011]. We generated ensembles of heterogeneous permeability distributions (x) according to a log normal distribution (x) . LN (µ, 2), as described in section 1.3.2. The statistical distribution is characterized by the geometric mean . = exp µ and variance 2. The investigation 3.3. Numerical Groundwater Flow Model layer L s s f range of . (l) (r) (h) km 0 - 45 2 600 2.0 0.30 100 - 1 000 100 100 1 000 ku 0 - 55 2 600 2.0 0.25 100 - 1 000 100 100 1 000 mo 80 2 800 2.0 0.20 5 - 200 5 20 200 mm 50 2 800 2.0 0.20 20 - 1 000 20 100 1 000 mu 100 2 800 2.0 0.20 5 - 200 5 20 200 so 100 2 400 2.7 0.05 0.02 - 2 0.02 0.2 2 sm 240 2 400 2.1 0.15 30 - 3 000 30 300 3 000 suB 160 2 400 2.1 0.13 30 - 300 30 100 300 suC 185 2 400 1.7 0.08 2 - 100 2 20 100 Table 12: Simulation parameters for sedimentary layers: mean thickness L [m], rock density s [kg m-3], thermal conductivity s [kgm2 s-3 K-1], porosity f [-], and permeability . [mD] (=[1e-15m2]). The ranges of . correspond to results found in the literature (Table 10.) (l), (r), and (h) are the three mean permeability values for every sedimentary layer, with abbreviations (l) for low, (r) for reference, and (h) for high. included three degrees of heterogeneity: mildly heterogeneous media, expressed by a variances of 2 =0.2, medium heterogeneous media with 2 = 1 and strongly heterogeneous media with 2 = 4. We used two different spatial correlation structure models. First, we applied a Gaussian covariance model, accordingtoEq.(1.18), characterizedby a horizontal and a vertical correlation length `h and `v. We generated heterogeneous fields with correlation lengths over four orders of magnitude: from local scale `h = 20 m (`v = 2 m) corresponding to the size of a grid cell, over regional scale `h = 200 m (`v = 20 m) and `h = 2000 m (`v = 50 m) up to basin scale `h = 20000 m (`v = 50 m) which corresponds to the maximal size of an undisturbed basin compartment. Second, we used the truncated power law covariance model, according to Eq.(1.19). In this model there is no finite correlation length, but correlation exists at any scale. We used an upper cut-off value of `max = 25000 m (`max = 50 m) which is above the hv maximal size of an undisturbed basin compartment in the Thuringian Basin. Every ensemble consists of 200 realizations with identical choice of statistical parameters. We performed simulations to investigate density dependent flow, where the fluid density can change due to dissolved salt or temperature variations. We applied a linear density model, according to Eq.(1.10): . = 0(1 + C C + T T ). We set the solute expansion coefficient to C =0.2, which refers to a maximal saturated brine density of C = 1200 kg/m3 for a normalized salt concentration C . [0, 1]. The thermal expansion coefficient is T = -0.0003 corresponding to minimal density of T = 991 kg/m3 for the maximal temperature of 38 C and a fresh water density of 0 = 1000 kg/m3 for 8 C (Table 11). We made some simplifying assumptions. We assumed the brine to be pure NaCl solution and neglected viscosity variations due to temperature and salt concentration. Furthermore, the simulation results cannot reproduce the real groundwater flow, which is naturally three dimensional. Nevertheless, the two dimensional approach allowed to investigate brine migration in shallow sedimentary basins and to gain insight into the role of geological fault structures, hydraulic parameters and fluid density variation on the flow, salt and temperature pattern. 3. Mechanisms of Salt Transport in the Thuringian Basin 3.4. Discussion of Simulation Results In the following we discuss, compare, and interpret the results for velocity distribution, salt concentration, and temperature pattern of different simulation scenarios.The results of the reference simulation serve as basis for comparison with all other simulation scenarios. We then present the results for changes in the mean permeability of different sedimentary layers and in the fault zones. It is followed by an interpretation of the simulation results with heterogeneous permeability distributions and an analysis of the impact of the statistical quantities. The final step involves the discussion of simulation results for density dependent flow. Besides visual comparison, we use quantitative characteristics to support interpretation results. ¯ For every simulation we computed the mean of the salt concentration C . [0, 1], as given in Eq.(1.12) which gives a relative value of the total amount of dissolved salt in in the entire domain. For every simulation scenario, we calculated the relative concentration difference C [%] and the relative temperature difference T [%] to the reference simulation, as given in Eq.(1.13). Solute and thermal Rayleigh numbers RaS and RaT are calculated according to Eq.(1.11) for density dependent flow. 3.4.1. Reference Simulation The initial point for all comparisons of different simulation scenarios is the reference simulation. Boundary and initial conditions as well as domain parameters are defined in section 3.3.3 and section 3.3.4. We assumed the permeability to be homogeneously distributed with mean values of all layers according to the reference values (r), defined in Table 12. The transport of salt and heat is conservative, meaning that a constant fluid density model is used, where salt and temperature do not impact on the flow pattern. Simulation results of velocity, normalized salt concentration distribution and the temperature distribution after 100000 years are visualized in Figure 15. We assume that these results represent the steady state solution because for all processes the changes to subsequent simulation time steps are less than 0.05%. The time needed to establish the regional flow pattern from the initial conditions can be neglected compared to the time needed to reach a steady state in temperature and salt concentration. The temperature distribution reaches steady state after 30000 years, whereas the salt requires approximately 100000 years to reach a stadium of stable salt stratification. The faster development of thermal steady state is due to the higher diffusivity of the heat transport process. The topography determines the flow pattern (Figure 15a), according to the choice of the constant hydraulic head flow boundary condition at the basin surface. Flow takes place from elevated regions to lower basin parts. The main area of recharge is the southwest basin margin at the Thuringian forest border fault. Discharge takes place locally, especially in flat basin regions and near the fault zones. We observe high velocities in the shallow Keuper layer (km/ku) and the Middle Buntsandstein (sm) aquifer due to the high permeabilities. We find very low velocity in the Röt aquiclude (so) and in the Lower Buntsandstein layers (suB/suC). Especially 3.4. Discussion of Simulation Results Figure 15: Simulation results for reference simulation: (a) Velocity distribution with arrows marking the flow direction, (b) normalized salt concentration with isolines at C = [0.1,...,0.9] and (c) Temperature distribution with isolines at distance of 2C. Sedimentary units are separated by bold black lines. Superelevation factor of 10. in the flat basin center flow is nearly stagnant, which corresponds to the findings of Meincke [1967]. The fault geometry has distinct effects on the flow pattern. Where the EichenbergGotha- Saalfelder fault zone and the Erfurter fault zone serve as transit areas for flow, the Sömmerdaer fault has a barrier effect on the flow causing the groundwater to rise. The salt concentration pattern (Figure 15b) corresponds to the distribution of velocities. The choice of permeabilities allow the development of large scale groundwater advection cell. They reach the deep aquifer units of the Buntsandstein and transport saline groundwater to nearsurface regions. Steep plumes of salt water develop in the two central undisturbed basin compartments and one in the vicinity of the Sömmerdaer fault. The mean concentration in the ¯¯ entire domain is C(k/m/s) =0.159, where the amount of salt in the shallow aquifer C(k/m) =0.016 is very low compared to the amount in the deep aquifer C¯(sm/su) =0.224 (Table 13). 3. Mechanisms of Salt Transport in the Thuringian Basin (k/m) (so) (sm/su) ¯ C(k/m/s) ¯ C(k/m) ¯ C(so) ¯ C(sm/su) C [%] r r r 0.157 0.014 0.030 0.224 0 h r r 0.154 0.004 0.025 0.221 2.66 l r r 0.171 0.022 0.028 0.239 3.73 r h r 0.133 0.010 0.013 0.190 6.16 r l r 0.194 0.012 0.050 0.274 7.88 r r h 0.189 0.023 0.031 0.266 7.82 r r l 0.142 0.009 0.022 0.202 5.22 r h h 0.163 0.014 0.013 0.231 4.44 r h l 0.124 0.007 0.012 0.177 6.44 r l h 0.194 0.016 0.050 0.273 8.73 r l l 0.179 0.008 0.046 0.254 4.70 h r h 0.172 0.006 0.023 0.247 6.56 h r l 0.140 0.002 0.019 0.203 5.16 l r h 0.174 0.028 0.027 0.257 7.51 l r l 0.143 0.014 0.023 0.202 5.26 h h r 0.133 0.005 0.012 0.191 6.05 l h r 0.139 0.013 0.015 0.197 5.92 h l r 0.175 0.002 0.033 0.266 7.35 l l r 0.218 0.032 0.075 0.299 10.03 h h h 0.149 0.007 0.013 0.215 3.31 h h l 0.122 0.002 0.001 0.177 6.44 h l h 0.172 0.002 0.026 0.248 7.16 h l l 0.178 0.001 0.037 0.256 5.22 l h h 0.165 0.019 0.016 0.233 4.55 l h l 0.126 0.011 0.013 0.180 6.36 l l h 0.216 0.044 0.073 0.292 10.80 l l l 0.178 0.018 0.054 0.262 5.33 ¯ Table 13: Mean concentration C, as defined in Eq.(1.12) and relative concentration difference C to the reference simulation, as defined in Eq.(1.13) for combinations of homogeneous permeabilities of the three hydraulic units shallow aquifer Keuper/Muschelkalk (k/m), aquiclude Röt (so) and deep aquifer Middle and Lower Buntsandstein (sm/su). The first three columns denote the choice of permeability : l=low, r=reference, h=high, according to the values given in Ta ¯ ¯¯ ble 12; C(k/m/s) is the mean concentration in the entire domain; C(k/m), C(so), and C¯(sm/su) are the mean concentrations in the corresponding hydraulic units. 3.4. Discussion of Simulation Results The distribution of temperature (Figure 15c) develops according to the thermal gradient of 30 Cper 1kmdepth: from 8 C at the surface to maximal 38.4 C at the bottom of the basin. The stable temperature stratification, which develops by conductive heat transport trough the solid and fluid phase, is disturbed by regional groundwater flow. Surface-near regions of lower and higher temperatures coincide with the areas of recharge and discharge: cold meteoric fresh water enters the system and is transported with the regional groundwater flow into deeper aquifer units. This can be observed especially in the environment of fault zones. In discharge areas the warmer groundwater from deeper regions is transported to the surface. This can be observed mainly in the flat undisturbed basin units, where also salt water plumes develop. Thus, conduction is the main transport mechanism in deep aquifer units, which is visible at the flat isolines. Whereas advection dominates the temperature distribution in the shallow aquifer units. 3.4.2. Scenarios of Homogeneous Permeability The combination of the three choices of permeability , as listed in Table 12, for every of the three hydraulic units results in 27 possible scenarios of mean permeability distribution. In Table 13 we give the calculated characteristic numbers of mean concentration C ¯ in the entire domain (k/m/s) and the hydraulic units Keuper/Muschelkalk (k/m), Röt (so) and Middle and Lower Buntsandstein (sm/su) as well as the relative difference to the reference simulation C for all simulations. To distinguish between the simulation scenarios, we use subscripts ¯ symbolizing the choice of . in the corresponding hydraulic unit; for instance Clrh describes the mean concentration for the scenario with low permeability in the shallow aquifer Keuper/ Muschelkalk (k/m)= l, reference permeability in the aquiclude Röt (so)= r and high permeability in the deep aquifer Middle/Lower Buntsandstein (sm/su)= h. The discussion and visualization of simulation results focus on the salt distribution of scenarios where permeabilities of only one hydraulic unit were modified compared to the reference simulation. For brevity, simulation results for the flow field are not visualized separately, since it can be inferred from the salt concentration distribution. Furthermore, the temperature distribution is not visualized and discussed, because changes to the reference temperature distribution are small. Variation of Permeability in the Shallow Aquifer The simulated concentration distribution for the scenarios with increased and decreased permeabilities in the Keuper and Muschelkalk layers are visualized in Figure 16. The overall flow pattern only changes little compared to the reference simulation. Number and location of salt plumes are similar, which is emphasized by the small concentration difference to the reference simulation of Chrr =2.66% and Clrr =3.73% (Table 13). We observe differences in the mean concentration of the entire domain which is higher for lower permeabilities and vice versa: C¯ lrr(k/m/s) = 0.171 and C¯ hrr(k/m/s) = 0.154 compared 3. Mechanisms of Salt Transport in the Thuringian Basin Figure 16: Distribution of normalized salt concentration with isolines at C = [0.1,...,0.9] of simulations with (a) increased and (b) decreased permeability in shallow aquifer Keuper/Muschelkalk (k/m). Sedimentary units are separated by bold black lines. Superelevation factor of 10. to C¯ rrr(k/m/s) =0.159. Significant differences occur in the shallow aquifer units Keuper and Muschelkalk (Table 13). The spreading of salt within the Muschelkalk is apparently damped for higher permeabilities. This effect can be explained by deeper intruding fresh water due to higher velocities. However, a change in permeability of the shallow aquifer has little impact on the flow pattern in the aquiclude Röt or the Buntsandstein aquifer. Variation of Permeability in the Aquiclude Simulated concentration distribution for varied permeabilities in the aquiclude Röt (so) are given in Figure 17. The changes compared to the reference simulation are significant, visible in Crhr =6.16% and Crlr =7.88% (Table 13). An increased mean permeability leads to more salt plumes with slightly changed position in comparison to those of the reference simulation. The salt plumes are steeper and have sharper interfaces. Furthermore, less salt is dissolved ( C¯ rhr(k/m/s) =0.133, Table 13). An increased permeability in the Röt allows a stronger communication between the shallow aquifer and the deep aquifer. Thus, the regional groundwater flow has a stronger impact in deeper basin units, resulting in more pronounced large scale advection cells. The flow and salt patterns change completely when reducing the permeability of the Röt. An impermeable Röt layer disconnects Muschelkalk and Middle/Lower Buntsandstein layers, 3.4. Discussion of Simulation Results Figure 17: Distribution of normalized salt concentration with isolines at C = [0.1,...,0.9] of simulations with (a) increased and (b) decreased permeability in aquiclude Röt (so). Sedimentary units are separated by bold black lines. Superelevation factor of 10. which results in a break down of large scale advection cells. In the flat basin center diffusive and dispersive salt transport dominates due to the very low groundwater flow velocities in this region. A strong salinization of the deep aquifer (sm/su) can be observed (Figure 17b). ¯ Although the mean salt concentration increased to Crlr(k/m/s) =0.194, the amount of salt reaching near-surface areas is reduced. Variation of Permeability in the Deep Aquifer Changes in the permeabilities of the deep aquifer units have strong impact on the flow pattern, visualized in in Figure 18. The concentration distribution for varied permeabilities in the deep aquifer differs significantly from the results of the reference simulation (Figure 15b). This is underlined by the high relative differences of Crrh =7.82% and Crrl =5.22% (Table 13). For higher permeabilities (sm/su)= h less salt plumes develop (Figure 18a). We observe only one large plume in the flat basin center. It shows that the occurrence of large scale groundwater advection is reduced. The salinization is stronger in the Lower Buntsandstein layers, visible in a higher amount of dissolve salt, C¯ rrh(k/m/s) =0.189 (Table 13). A reduction of permeability in the Buntsandstein aquifer has contrary effects. The concentration distribution shows more and smaller salt plumes than the reference simulation (Figure 18b compared to Figure 15b). The total amount of dissolved salt is reduced, C¯ rrl(k/m/s) =0.142, as well as the mean concentration for all hydrological units (Table 13). 3. Mechanisms of Salt Transport in the Thuringian Basin Figure 18: Distribution of normalized salt concentration with isolines at C = [0.1,...,0.9] of simulations with (a) increased and (b) decreased permeability in deep aquifer Middle and Lower Buntsandstein (sm/su). Sedimentary units are separated by bold black lines. Superelevation factor of 10. Variation of Permeabilities in the Faults The impact of the permeability distribution in the fault zones is investigated by multiplying the mean permeabilities of all sedimentary layers in the fault zones with a certain factor, while keeping the mean permeability distribution in the undisturbed basin compartments constant. Simulation results are visualized in Figures 19 and 20. Table 14 contains the calculated mean ¯ concentration C in the entire domain and in the sedimentary units, as well as the relative difference to the reference simulation C. The changes in flow and salt patterns are recognizable but overall small. The number of salt plumes and the total amount of dissolved salt is constant. This is underlined by the small concentration differences to the reference simulation C for all scenarios, which are between 1.4% and 4% (Table 14). We observe differences to the reference simulation locally in the shape and positions of the salt plumes. Most significant are the changes in flow pattern for highly impermeable faults, Figure 20b. Differences can be seen in the location of the salt plume at the Sömmerdaer fault and of the plume in the vicinity of the Eichenberg-Gotha-Saalfeld fault zone. Furthermore, the salt distribution in the shallow aquifer of the central salt plume is modified. These changes can be explained by differences in the velocity of inflowing fresh water in the fault. 3.4. Discussion of Simulation Results Figure 19: Distribution of normalized salt concentration with isolines at C = [0.1,...,0.9] of simulations with increased permeability . in fault zones: (a) 10 · , (b) 100 · . Sedimentary units are separated by bold black lines. Superelevation factor of 10. Figure 20: Distribution of normalized salt concentration with isolines at C = [0.1,...,0.9] of simulations with decreased permeability . in fault zones: (a) 0.1 · , (b) 0.01 · . Sedimentary units are separated by bold black lines. Superelevation factor of 10. 3. Mechanisms of Salt Transport in the Thuringian Basin * ¯ C(k/m/s) ¯ C(k/m) ¯ C(so) ¯ C(sm/su) C [%] 1 0.159 0.016 0.030 0.224 0 10 0.158 0.016 0.033 0.222 1.43 100 0.154 0.016 0.031 0.217 1.65 0.1 0.159 0.013 0.029 0.225 1.97 0.01 0.156 0.013 0.025 0.221 4.04 ¯ Table 14: Mean concentration C, as defined in Eq.(1.12) and relative concentration difference C to the reference simulation, as defined in Eq.(1.13) for different permeabilities in fault zones: * denotes the factor of modified permeability in fault zones, ¯ as defined in section 3.3.4; C(k/m/s) is the mean concentration in the entire domain, C¯(k/m), C¯(so) and C¯(sm/su) are the mean concentrations in the corresponding hydraulic units. 2 `h [m] `v [m] ( ¯ C(k/m/s)) ( ¯ C(k/m)) ( ¯ C(so)) ( ¯ C(sm/su)) hC) hC) 0 0 0 0 0.159 0.014 0.030 0.224 0 0 10 0.2 20 2 0.153 0.014 0.028 0.217 1.15 0.84 20 0.2 200 20 0.156 0.014 0.030 0.220 2.27 1.39 30 0.2 2 000 50 0.157 0.014 0.028 0.222 3.59 1.96 40 0.2 20 000 50 0.159 0.014 0.027 0.225 3.63 2.26 50 0.2 25 000 50 0.158 0.014 0.027 0.223 3.08 1.86 11 1 20 2 0.156 0.017 0.031 0.219 2.39 1.77 21 1 200 20 0.165 0.020 0.037 0.231 4.79 3.45 31 1 2 000 50 0.170 0.020 0.036 0.238 6.73 4.26 41 1 20 000 50 0.170 0.017 0.032 0.240 6.32 4.20 51 1 25 000 50 0.168 0.017 0.031 0.236 5.52 3.61 12 4 20 2 0.158 0.022 0.038 0.220 5.11 3.93 22 4 200 20 0.169 0.026 0.049 0.232 7.41 5.64 32 4 2 000 50 0.188 0.030 0.058 0.259 10.49 7.39 42 4 20 000 50 0.190 0.028 0.054 0.262 9.91 6.95 52 4 25 000 50 0.181 0.025 0.048 0.251 8.84 6.20 Table 15: Ensemble mean concentrations and relative concentration differences for all ensembles of heterogeneous permeability distributions: The first column denotes the ID of the ensemble. The choice of statistical parameters is given by the variance 2 and the horizontal and vertical correlation lengths `h and `v, respectively. hC¯(k/m/s)) is the ensemble mean concentration in the entire domain, hC¯(k/m)i, hC¯(so)i, and hC¯(sm/su)) are the ensemble mean concentrations in the corresponding hydraulic units. They are calculated as ensemble average over the mean concentrations C ¯ of the 200 realization; hC) [%] is the ensemble average of the relative differences of concentration Ci, i =1,..., 200 between single realization i and the reference simulation. Whereas hC) [%] is the relative difference of concentration between the ensemble mean hC) and the reference simulation. 3.4. Discussion of Simulation Results 3.4.3. Scenarios of Heterogeneous Permeability Distribution For every ensemble, we calculated the spatially distributed ensemble concentration hC) as the mean over the concentration distribution of all realizations. Table 15 contains the calculated quantitative measures for the simulated ensembles of heterogeneous permeability distributions. The ensemble mean concentration hC¯) is the ensemble average of the mean concentrations C ¯ of all realizations. The relative ensemble difference hC) is the ensemble average of the relative difference C between every single realization and the reference simulation, whereas the ensemble relative difference in concentration hC) denotes the difference between the ensemble mean hC) and the reference simulation. Impact of the Correlation Structure We compare the simulation results for different correlation structures while keeping the variance at a medium level of 2 = 1. First, we investigate the impact of a Gaussian correlation structure for four cases of correlation length (`h = 20 m, 200m, 2000 m and 20000 m) ranging from local to basin scale, as defined in section 3.3.4. Second, we analyze an ensemble with a truncated powerlawcovariancestructurewithcutoff alengthof 25000 m. Figure 21 shows the simulated salt concentration distributions for two realizations and the ensemble mean of ensemble 11. The correlation length of `h = 20 m ranges in the size of one grid cell. It can be seen that a small correlation length causes local fluctuations in the salt concentration, but does not impact on the shape and the position of the salt plumes compared to the homogeneous reference simulation (Figure 15b). This is underlined by the small mean difference of hC11) =2.39% (Table 15). The differences between realizations (Figure 21a and 21b) are little. Thus, the ensemble mean concentration pattern hC11) is nearly identical to the concentration distribution with homogeneous permeabilities (Figure 21c compared to Figure 15b). This fact is underlined by the small relative difference hC11) =1.52% and a similar amount of dissolved salt hC¯ 11(k/m/s) ) =0.156 compared to the reference simulation. We visualize the simulation results for a medium scale correlation length of `h = 200 m (ensemble 12) in Figure 22. Small changes exist in position and shape of salt plumes for the heterogeneous realizations (Figures 22a and b) compared to the homogeneous reference simulation (Figure 15b). Local fluctuations are more pronounced in comparison to the results for a small correlation length of 20 m. This is supported by a higher relative concentration difference of hC21) =4.79%. An increased spreading of salt is visible in a higher amount of dissolved salt hC¯ 21(k/m/s) ) =0.165 (Table 15) compared to the reference simulation. The ensemble mean of heterogeneous simulations hC21) in Figure 22c mirrors changes in the concentration pattern: salt plumes are damped and there exists a certain range in which salt plumes develop. 3. Mechanisms of Salt Transport in the Thuringian Basin Figure 21: Distribution of normalized salt concentration with isolines at C = [0.1,...,0.9] of simulations with heterogeneous permeability distribution: (a) and (b) are single realizations and (c) is the mean over 200 realizations of ensemble 11 with small correlation length (. = 20 m). Sedimentary units are separated by bold black lines. Superelevation factor of 10. 3.4. Discussion of Simulation Results Figure 22: Distribution of normalized salt concentration with isolines at C = [0.1,...,0.9] of simulations with heterogeneous permeability distribution: (a) and (b) are single realizations and (c) is the mean over 200 realizations of ensemble 21 with medium correlation length (. = 200 m). Sedimentary units are separated by bold black lines. Superelevation factor of 10. 3. Mechanisms of Salt Transport in the Thuringian Basin Figure 23: Distribution of normalized salt concentration with isolines at C = [0.1,...,0.9] of simulations with heterogeneous permeability distribution: (a) and (b) are single realizations and (c) is the mean over 200 realizations of ensemble 31 with large correlation length (. = 2000 m). Sedimentary units are separated by bold black lines. Superelevation factor of 10. Figure 23 shows the simulation results for large scale correlation length of `h = 2000 m. Local fluctuation are of larger scale and can lead to significantly different salt patterns for different realizations (Figure 23a and b). We observe differences in number, position and shape of salt plumes. The differences to the reference simulation (Figure 15b) are large, supported by the relative ensemble difference of hC31) =6.73%. A higher amount of dissolved salt hC¯ 31(k/m/s)) =0.17 indicates a wider spread of salt plumes. The ensemble concentration pattern hC31) (Figure 23c) shows a flat salt distribution without clear salt plumes. It displays the range of positions of salt plumes. In Figure 24, we visualize the simulation results for a correlation length of `h = 20000 m, which corresponds to the maximum distance of an undisturbed basin compartment. The results show significant differences in number, position, and shape of salt plumes. Realizations differ strongly 3.4. Discussion of Simulation Results Figure 24: Distribution of normalized salt concentration with isolines at C = [0.1,...,0.9] of simulations with heterogeneous permeability distribution: (a) and (b) are single realizations and (c) is the mean over 200 realizations of ensemble 41 with correlation length at basin scale (. = 20000 m). Sedimentary units are separated by bold black lines. Superelevation factor of 10. among each other and from the homogeneous reference simulation, visible at hC41) =6.32% (Table 15). Due to the very large correlation length, the results can be similar to salt patterns of simulations with homogeneous permeabilities, as discussed in section 3.4.2. The amount of dissolved salt is higher (hC¯ 41(k/m/s)) =0.17) compared to homogeneous simulations. The pattern of the ensemble mean hC41) (Figure 24c) is similar to that of the large scale correlation length `h = 2000 m (Figure 23c) mirroring the areas of salt plumes. Figure 25 shows the simulation results with the truncated power law correlation model, displaying a heterogeneous medium with correlation on all scales. In the resulting concentration patterns, we observe local fluctuations, as in simulations with small correlation lengths and changes in position and shape of salt plumes, as observed in simulations with high and basin scale correlation lengths. The concentration distributions differ in average hC51) =5.52%. 3. Mechanisms of Salt Transport in the Thuringian Basin Figure 25: Distribution of normalized salt concentration with isolines at C = [0.1,...,0.9] of simulations with heterogeneous permeability distribution: (a) and (b) are single realizations and (c) is the mean over 200 realizations of ensemble 51 with infinite correlation length (truncated power law model). Sedimentary units are separated by bold black lines. Superelevation factor of 10. The mean concentration pattern hC91) (Figure 25c) mirrors changes compared to the homogeneous reference simulation, as for the large scale correlation length. Generally, the resulting salt concentration patterns are very similar to those with correlation length f = 2000 m and `= 20000 m, which shows that the effect of large correlation lengths is dominant. Impact of the Variance We analyze simulations with three choices of variance: mildly heterogeneous media (2 = 0.2), medium heterogeneous media (2 =1), and strongly heterogeneous media (2 =4). Figure 26 shows the simulated salt distributions for three realization with different variances and a Gaussian correlation length of `h = 2000 m. We chose the spatial correlation structure 3.4. Discussion of Simulation Results Figure 26: Distribution of normalized salt concentration with isolines at C = [0.1,...,0.9] of simulations with heterogeneous permeability distribution for different variances: realizations from ensembles 30, 31 and 32 with identical correlation structure (`h = 2000 m) and (a) 2 =0.2, (b) 2 =1 and (c) 2 =4. Sedimentary units are separated by bold black lines. Superelevation factor of 10. identically for all variances, which allows a direct comparison of realizations. It can be seen that a higher variance goes along with a stronger salinization. This is valid for all choices of correlation structure, which is visible in increasing values of hC¯) for increasing variance for all correlation lengths (Table 15). However, the variance does not directly impact on number and position of salt plumes, only on the shape due to the higher amount of dissolved salt. This is in particular valid in combination with larger correlation lengths. Heterogeneity contributes to the mechanical dispersion [Bear, 1972]. The variance controls the range of permeability values, amplifying differences between areas of higher and lower permeability. As a consequence, the larger the variance, the stronger are the variations in velocity. Higher velocities result in increased mechanical dispersion and thus in higher amounts of dissolved salt. This effect can be observed in particular in the shallow aquifer units. 3. Mechanisms of Salt Transport in the Thuringian Basin density model ¯ C(k/m/s) ¯ C(k/m) ¯ C(so) ¯ C(sm/su) C [%] T [%] 0 0.157 0.014 0.030 0.224 0 0 (T ) 0.157 0.014 0.029 0.222 1.00 0.17 (C) 0.309 0.079 0.149 0.408 15.87 2.17 (T, C) 0.308 0.078 0.146 0.406 15.71 2.14 Table 16: Mean concentration C¯, as defined in Eq.(1.12), relative difference in concentration C and relative difference in temperature T to the reference simulation, as defined in Eq.(1.13) for different fluid density models: 0 is the constant density model, (T ) is the purely temperature dependent density model, (C) is the purely salt concentration dependent density model, and (T,C) is the concentration and temperature dependent density model according to Eq. (1.10). C¯(k/m/s) is the mean concentration in the entire domain, C¯(k/m), C¯(so) and C¯(sm/su) are the mean concentrations in the corresponding hydraulic units. 3.4.4. Density Dependent Flow Impact of the Temperature The simulation results for a temperature dependent fluid density model (T ) are given in Figure 27. The salt and temperature distributions resemble those of the constant density reference simulation in Figure 15. We observe no differences in the temperature distribution visually. Small changes in the salt pattern can be found in the form of the central salt plume. The marginal differences are underlined by the nearly identical values of the mean concentration (Table 16) and the small values of the relative difference in temperature T =0.17% and salt concentration C =1.0%. The marginal changes for a temperature dependent density model are due to the fact that the temperature differences are too small to initiate convective flow. According to the linear density model described in Eq.(1.10), the density of water can vary between 1000 kg/m3 and 991kg/m3, which corresponds to the thermal expansion coefficient of T = -0.0003 and the minimal and maximal temperature of 8 C and 38.4 C respectively. The maximal thermal Rayleigh number RaT, as defined in Eq.(1.11) for all sedimentary units can be found in the Middle Buntsandstein aquifer with RaT = 11.7. This is beneath the critical value of 4p which characterizes the onset of convection [Nield and Bejan, 1999]. Even for higher permeabilities, thus higher thermal Rayleigh numbers, regional groundwater flow would superpose convective processes due to higher flow velocities. From the Rayleigh theory we know that temperature induced convection requires highly permeable media of sufficiently large vertical extent and which is undisturbed from other flow sources. Since this is not present in the Thuringian Basin, conduction and advection dominate the heat transport process. Thus, temperature dependent density driven flow is negligible in the Thuringian Basin due to the significant regional flow. 3.4. Discussion of Simulation Results Figure 27: Simulated temperature and salt distribution for flow impacted by temperature induced density variations and regional groundwater flow: (a) salt concentration distribution with isolines at C = [0.1,...,0.9] and (b) temperature distribution with isolines at distance of 2 C. Sedimentary units are separated by bold black lines. Superelevation factor of 10. 3. Mechanisms of Salt Transport in the Thuringian Basin Figure 28: Results of temperature and salt distribution for flow impacted by density variations due to salt concentration differences and regional groundwater flow: (a) salt concentration distribution with isolines at C = [0.1,...,0.9] and (b) temperature distribution with isolines at distance of 2 C. Sedimentary units are separated by bold black lines. Superelevation factor of 10. Impact of Salt Concentration Effects of dissolved salt on the fluid density result in salt concentration and temperature patterns given in Figure 28. Significant differences are visually recognizable compared to the reference simulation in Figure 15b. More salt is dissolved for a salt dependent density model (C) than for the constant density model 0: C¯(k/m/s) =0.309 compared to C¯(k/m/s) =0.159, respectively (Table 16). The salt plume pattern is less pronounced. High concentration brine (C> 0.7) does not move upwards. Salt is stronger laterally distributed, in particular in the central basin compartments. Changes in temperature distribution (Figure 28b) are less pronounced but observable compared to the reference simulation (Figure 15c). The temperature isolines are smoother, which results from changed velocity fields in the deep aquifer, due to salt mass effects. The relative difference in temperature pattern amounts to T =2.17%. The strong impact of salt on the flow and salt pattern can be explained by the interplay of density differences and the regional groundwater flow. Full salt saturation can lead to a fluid density of up to 1200 kg/m3. This results in large solute Rayleigh numbers RaS in all sedimentary layers. In the deep aquifer the values range between RaS = 800 - 6000, where the highest values occur in the Middle Buntsandstein of the undisturbed central basin compartment. In areas of discharge water pressure gradients cause an upward movement of 3.4. Discussion of Simulation Results Figure 29: Results of temperature and salt distribution for flow impacted by density variations due to temperature and salt concentration differences and regional groundwater flow: (a) salt concentration distribution with isolines at C = [0.1,...,0.9] and (b) temperature distribution with isolines at distance of 2 C. Sedimentary units are separated by bold black lines. Superelevation factor of 10. brine. The effect of an increased fluid density due to solved salt leads to the counteracting effect of downward convection. It results in smoothed pressure gradients in horizontal direction but increased vertical velocities. Higher velocities cause an increase in dispersion. Furthermore, local eddies develop as described in Holzbecher et al. [2010], which enhance local mixing. Both effects result in a stronger spreading of salt. Impact of Temperature and Salt Concentration Simulation results for a salt and temperature dependent density model (C,T) are given in Figure 29. The salt and temperature distributions resemble the patterns with a purely salt dependent density (C) (Figure 28). Changes in the characteristic numbers are small: the total amount of salt and the relative difference to the reference simulation are marginally lower than for purely salt dependent density (Table 16). In general, the effect of salt is slightly reduced, but still dominant. The reduction can be explained by the opposite character of salt and of temperature on fluid density. 3. Mechanisms of Salt Transport in the Thuringian Basin 3.5. Conclusions The aim of the study was to investigate the main fluid dynamics within the Thuringian Basin. We performed transient coupled simulations of hydraulic, thermal, and mass transport processes on a selected profile to quantify the impact of the permeability distribution and density differences on the flow pattern, the salt concentration and the temperature distribution. To capture the uncertainty in permeability, which determines regional groundwater flow, we performed simulations with different mean permeabilities in the main hydraulic units as well as in the fault zones. We further analyzed effects of heterogeneity. Thereby permeability was characterized by three parameters: mean, variance and correlation structure. In a first step, we examined the effect of mean permeability varying within a realistic range by conducting simulation with different distributions of homogeneous permeability within the three hydraulic units shallow aquifer (Keuper/Muschelkalk), aquiclude (Röt), and deep aquifer (Middle and Lower Buntsandstein). In a second step, we investigated the impact of heterogeneity using ensembles of heterogeneous permeability distributions with different values of variance and correlation structure. Thirdly, we examined the impact of salt mass and temperature on the groundwater flow. We performed simulations with a non-constant fluid density model, depending on salt concentration and temperature differences. The main driving force of fluid flow in the Thuringian Basin is the water pressure gradient, due to differences in the topography. The shallow basin structure causes the regional groundwater flow to impact on flow directions and velocities down to the deep aquifer units, connecting to the Zechstein salt in up to 1000 m depth. Especially in vicinity of faults, meteoric fresh water can infiltrate deeply into the basin. It results in high flow velocities and causes anomalies in the temperature stratification. Areas of discharge are the flat undisturbed basin compartments located in the basin center. Since flow velocities are proportional to the permeabilities of the sedimentary layers, their distribution significantly determines flow and salt concentration pattern. Simulation results with permeability variations in all hydraulic units indicate that variations in the permeability of the shallow aquifer (Keuper/Muschelkalk) have little impact on the salt distribution in the deep aquifer (Middle and Lower Buntsandstein). On the contrary, the permeability of the aquiclude (Röt) determines the amount of communication between deep and shallow aquifers and thus the flow pattern: the lower the permeability, the lower is the interaction and the more salt is dissolved in the deep aquifer, but the less salt is transported to surface near regions. Changes in permeability of the deep aquifer have significant impact on the amount of dissolved salt: the higher the permeability, the more salt is dissolve, but the salt plume pattern is less pronounced, thus less salt is transported to the surface. Heterogeneity in permeability can have significant impact on the developing flow pattern and thus the salt distribution. In general, an increase in correlation length and variance goes along with a deviation of the flow and salt pattern from results for a homogeneous permeability distribution. Fluctuations in the salt concentration pattern correspond to the choice of correlation length: small correlation lengths in the range of 20m cause local fluctuations with similar salt 3.5. Conclusions patternsashomogeneoussimulation. Largecorrelationlengthsintherangeof 2000 mcontrol large ranges of hydraulic conductivity. Thus heterogeneous simulation results can vary significantly in number and position of salt plumes, as it was the case for homogenous simulations with increased and decreased conductivity in the entire sedimentary layers. The truncated power law correlation structure could capture the effects of correlation lengths on all scales. A comparison with the Gaussian correlation structure for different correlation lengths showed that effects of large scale correlation lengths dominate over those of small scale correlation lengths. The variance as degree of heterogeneity has a different effect. Differences in permeability are stronger the larger the variance is. The resulting spread in local flow velocity enhances mechanical dispersion. Thus, higher variances only slightly change the overall flow and salt pattern but strongly impact on the amount of dissolved salt. The effect of an increased amount of dissolved salt is also given for increasing correlation length, although it is not that strong as for the variance. In cases of large correlation lengths and large variances, the deviation from salt distributions for homogeneous permeabilities is the strongest. Effects of a non-constant fluid density can impact on the flow pattern and thus the distribution of salt and temperature significantly. This is in particular the case for density variations due to salt concentration differences. High amounts of dissolved salt result in an increased density which causes downward movement of fluid. The counteracting effect of regional groundwater flow in discharge areas changes the velocity field and leads to the development of local eddies. It results in an enhanced mixing and thus an increase in the amount of dissolved salt. Solute Rayleigh numbers above the critical value of 4p in all aquifer units confirm these findings. On the contrary, the impact of temperature is little. Thermal induced convection is not present due to low temperature gradients and high regional flow velocities. These results from simulations are supported by Rayleigh theory: calculated thermal Rayleigh number RaT are lower than the critical value for onset of convection. The numerical study showed that in general two main transport mechanism interact in the Thuringian Basin: topography driven regional groundwater flow and convection due to fluid density differences. The interaction of both processes results in a wide spread of salt in the deep aquifer (Middle and Lower Buntsandstein) in the central basin. Salt transport to near-surface regions can be observed locally, especially in the basin center. In contrast to the dynamics of deep sedimentary basins like the North East German Basin, the shallow character of the Thuringian Basin inhibits the evolution of large thermal convection cells. From Rayleigh theory it is known that onset of convection requires thick, high permeable aquifers which are undisturbed from external flow. Those units are not present in the Thuringian Basin. However, other mechanism are responsible for brine transport to near-surface regions, which we showed by this numerical study. 3. Mechanisms of Salt Transport in the Thuringian Basin Part III. Determining Aquifer Heterogeneity on Local Scale 88 4. The Extended Thiem’s Solution -Including the Impact of Heterogeneity 4.1. Introduction Determining hydraulic properties of an aquifer has been a matter of research for several decades. Not only for characterizing groundwater flow but also for describing transport processes in the subsurface a good perception of the heterogeneous structure of porous media is necessary, see for e.g. Dagan [1989]; Gelhar [1993]; Rubin [2003]. Pumping tests are a widely used tool to identify hydraulic parameters which effect the groundwater flow pattern. Analyzing Darcy’s Law in combination with the continuity equation for steady state pumping tests in homogeneous porous media leads to the well known Thiem’s solution hThiem(r)= - Qw ln r + h(R) . (4.1) 2LK R It describes the drawdown of the hydraulic head h(r) depending on the radial distance from the well r for homogeneous hydraulic conductivity K. Thiem’s solution is valid in a confined aquifer of thickness L with fully penetrating well and the constant discharge Qw; h(R) is a known reference head at the distance R from the well. The applicability of Thiem’s solution to pumping tests in heterogeneous media is limited. It requires a representative conductivity value K for the whole range of the depression cone. As stated by Matheron [1967], a single representative K-value for well flow does not exist, due to the emergence of different mean conductivities characterizing the behavior near and far from the well. Since then an enormous amount of work has been devoted to find representative conductivity descriptions for well flow and to estimate statistical parameters of the hydraulic conductivity from drawdown data. For a detailed review see Sánchez-Vila et al. [2006]. Most of the studies are limited to a two dimensional analysis, like Shvidler [1966]; Desbarats [1992]; Sánchez-Vila et al. [1999]; Copty and Findikakis [2004]; Neuman et al. [2004, 2007]; Dagan and Lessoff [2007]; Schneider and Attinger [2008] and many more. Only a few authors addressed the impact of modeling the conductivity in three dimensions upon radial flow [Indelman and Abramovich, 1994; Indelman et al., 1996; Guadagnini et al., 2003]. In particular, only Firmani et al. [2006] presented a fully three dimensional numerical investigation. In order to find a description of the hydraulic head field in pumping tests for heterogeneous media, the conductivity K(x) is commonly modeled as a log-normal distributed spatial random function. Based on this assumption Indelman and Abramovich [1994] solved an averaged Darcy’s Law and presented a fundamental solution for the mean head distribution for arbitrary boundary conditions. Since it is given in Fourier space, only approximate solutions in real space are available. In this line Indelman et al. [1996] performed a perturbation expansion in the variance 2 of ln K(x) to present a first order solution in the hydraulic head for well flow. The result has been expanded by several authors to higher orders, e.g. Fiori et al. [1998]; Indelman [2001]. Additionally, Guadagnini et al. [2003] presented a three-dimensional steady state 4. The Extended Thiem’s Solution -Including the Impact of Heterogeneity solution for mean flow based on recursive approximations of exact non-local moment equations. The implicit character of these head solutions inhibits the application to analyze pumping test drawdowns directly. Furthermore, the reliability of a perturbation approach to describe well flow is questionable due to a ergodicity breakdown near the well which is mathematically a singularity. None of the solutions could reproduce the head drawdown of a pumping test exactly as numerical investigation showed [Guadagnini et al., 2003; Firmani et al., 2006] nor allowed an inverse estimation of the parameters of K(x) in highly heterogeneous porous media. Making use of the equivalent conductivity Keq as defined by Matheron [1967] and their first order solution, Indelman et al. [1996] derived the expression Keq(r)= Kwell(1-(r))+(r)Kefu. It relates the near well representative conductivity Kwell to the far field value Kefu (in detail discussed in section 4.2.2) by a weighting factor (r) which depends on the statistical parameters of K(x) and the radius r. This description was used by Firmani et al. [2006] for inverse parameter estimation from numerical pumping tests. But as their simulations showed, the description of Keq(r) is only valid for small variances 2 up to 0.5. Furthermore, the estimation of the parameters is of high uncertainty. To overcome the above mentioned limitations Schneider and Attinger [2008] showed that in two dimensions another description for the conductivity, respectively transmissivity, is appropriate to describe well flow effectively. They introduced a new approach by applying an upscaling technique to the flow equation to derive their representative description for the transmissivity T CG(r), depending on the radial distance and the statistical parameters. Based on that they performed forward simulations to achieve a head drawdown from T CG(r) and compared it to ensemble averages of simulated two dimensional pumping tests in heterogeneous media. They stated that this method allows a much better parameter estimation for T (x) than existing methods do. In this study we do not only extend the results of Schneider and Attinger [2008] to three dimensions but will go one step further by introducing a closed form solution for the effective well flow hydraulic head hefw(r). It describes the depression cone of a three dimensional pumping test in heterogeneous media effectively. This new solution hefw(r) can be understood as an extension of Thiem’s formula (4.1) to heterogeneous media. It accounts for the statistical parameters of K(x) and reproduces the vertical mean hydraulic head field at every radial distance r from the well preserving the flow rates. In contrast to existing solutions, hefw(r) does not result from a perturbation analysis of the mean head by expansion on the variance 2. Therefore, it is also valid for highly heterogeneous media. Furthermore, hefw(r) directly allows to estimate parameters of K(x) without the detour to a representative description of the conductivity as done by Indelman et al. [1996] and Firmani et al. [2006]. After stating the problem, we will shortly sum up known results for near and far field representative conductivities of pumping tests in section 4.2. In part 4.3, we introduce the upscaling method Coarse Graining and apply it to three dimensional well flow resulting in the representative conductivity description KCG(r). In section 4.4, we derive hefw(r) by solving the radial 4.2. Statement of the Problem flow equation with KCG(r) and perform a sensitivity analysis for the parameters of K(x) on the drawdown hefw(r). We finally prove the applicability of hefw(r) by analyzing three dimensional numerical pumping tests in highly heterogeneous media. Moreover, we implement an inverse estimation procedure to infer on the statistics of K(x) in part 4.5. The results presented within this section refer to the publication Zech et al. [2012]. 4.2. Statement of the Problem 4.2.1. Flow and Conductivity Model In this study we focus on steady state pumping tests with fully penetrating well in a confined aquifer, where water is extracted at a constant pumping rate Qw. The three dimensional flow is characterized by Darcy’s Law in combination with the continuity equation -r(K(x)rh(x)) = Q(x) , (4.2) with x =(x1,x2,x3). The sink/source term Q(x) can be written as Qw(xw - x), using the delta distribution function , where xw is the position of the well. To cover the spatial structure of the hydraulic conductivity K(x) we model it as log-normal distributed field K(x) . LN (µ, 2), meaning that ln K(x)/K0 is a normal distributed quantity with mean µ and variance 2. The spatial correlation structure described by a Gaussian shaped covariance model CV(x - x0)= 2(x - x0) includes the correlation function ! 222 x1 x2 x3 (x) = exp-- - , (4.3) `2`2`2 123 where `i is the correlation length in ith direction. To reflect the anisotropic structure of three dimensional heterogeneous media, we presume the correlation length in both horizontal directions to be identical `1 = `2 = g whereas in vertical direction we assume a smaller correlation length `3 = e`, with e . [0, 1] denoting the anisotropy ratio. 4.2.2. Far and Near Field Conductivities The issue of representative conductivity values for well flow has been discussed intensively over several years, starting with Shvidler [1966] and Matheron [1967]. However, most of the studies are limited to a two dimensional analysis. For three dimensional convergent flow Indelman and Abramovich [1994] concluded that the far field behavior is best covered by the effective hydraulic conductivity for uniform flow in three dimensions 1 Kefu = KG exp 22 - (e) , (4.4) 4. The Extended Thiem’s Solution -Including the Impact of Heterogeneity Figure 30: Relationship between Kefu, KA and KH. The values on the left refer to a variance of 2 =1 and KG =1m/s and are drawn to scale. where (e) is the anisotropy function, known from Dagan [1989], e 1 q (e)= parctan( 1/e2 - 1) - e. (4.5) 2(1 - e2)1 - e2 Depending on the degree of anisotropy, (e) varies between  13 (e=1) and zero (e=0), thus causing Kiso efu = KG exp  16 2to be the limit for isotropic media and the arithmetic mean KA =  2to be the limit for stratified media. KG exp 12 For a representative description of the conductivity near the well, in the following denoted by Kwell, different results emerge in literature (e.g. Indelman et al. [1996]; Indelman and Dagan [2004]), depending on the description of the discharge at the well. Theoretically, either a constant head hw (corresponding to a Dirichlet boundary condition) or a constant pumping rate Qw (Neumann boundary condition) can be applied. A constant pumping rate, in the following denoted as BCH, leads to a varying head along the well resulting is the harmonic mean of the hydraulic conductivity values at the well. On the contrary, a constant head gives the arithmetic mean as representative near well conductivity value. If ergodicity is fulfilled the  means are given by KA = KG exp 12 2 and KH = KG exp- 12 2 . An additional approach was introduced by Indelman et al. [1996] where a constant discharge was subdivided into fluxes proportional to the local conductivities along the well. This assumption leads to similar conditions as the constant head boundary condition, giving Kwell = KA, we will therefore call this BCA. All of these results were confirmed by numerical simulations by Firmani et al. [2006]. An important feature of three dimensional well flow is the asymmetric relation between KA, KH and Kefu, especially when taking anisotropy into account. As shown in Figure 30 the distance between Kefu and KA is much smaller than between Kefu and KH, even becoming zero for stratified media. Unlike in two dimension where we have K2d = KG and thus efu 23 | ln(K2d efu/KH)| =efu/KA)| = | ln(K2d 3d(iso) whereas | ln(K/KH)| efu = 2 . 12 2, we find in three dimensions | ln(KThis fact becomes important when comparing results for 3d(iso) efu /KA)| = 2 the different boundary conditions BCA and BCH at the well. 4.3. Method of Coarse Graining 4.3. Method of Coarse Graining We use the upscaling technique Coarse Graining as introduced by Attinger [2003] to gain a representative description for the well flow conductivity. Schneider and Attinger [2008] already applied Coarse Graining to radial convergent flow in two dimensions, focusing on regional scale flow. The following procedure of deriving KCG(r) will be similar, but we focus on three dimensional media, additionally incorporating anisotropy. 4.3.1. Coarse Graining as Spatial Filtering Procedure The idea of the Coarse Graining method is to apply a spatial filter of variable volume size on the flow equation (4.2) to transform it to a coarser scale. The approach was originally developed for Large Eddy Simulations, see Layton [2002]. The background of using Coarse Graining in hydrogeological modeling is to find a representative conductivity on a coarser scale still considering sub-scale effects. The procedure is based on averaging the flow equation over a filter volume controlled by the filter length scale . Mathematically this is expressed by a convolution with a smoothing function f(x). It results in a flow equation on coarser scale for the filtered heads h(x)= R 1 R 0R3 f(x0)h(x' + x)d3x' where the effects of head fluctuations smaller than . are f(x0)dx covered by the scale dependent -filtered conductivity KCG(x). Recapitulating the method . derived in Attinger [2003], the head equation is transformed to Fourier space, then the filtering procedure is applied to the transformed head equation. After some mathematical treatment and certain approximations the equation is transformed back and we result in the filtered head  KCG -r(x)rh(x)= Q(x) , (4.6) . where Q(x) is the filtered sink term. KCG(x)= KCG()+ K˜ CG(x,) denotes the -filtered . A hydraulic conductivity, which is still a spatial random function on a coarser scale. It is composed of the filtered mean conductivity KCG(), depending on , but not on x and the filtered A fluctuation part K˜ CG(x,). The result for KCG() depends on the geostatistical model and on the choice of the filter A  x function f(x). Using a Gaussian shaped filter function f(x) /rexp22 and making use of renormalization theory, Attinger [2003] presented a closed form solution for isotropic media  by KCG()= Kiso 321+ 2 - 3 efu is defined in Eq.(4.4) with e = 1. A efu exp 14`22 , where Kiso 4.3.2. Coarse Graining for Well Flow Coarse Graining comprises the possibility of applying a filter of variable volume size on the flow equation which is of significant importance when dealing with well flow. For pumping tests the singular character of the source causes a strong influence of local heterogeneities on the flow pattern near the well, where the impact decreases with distance from the well. We therefore 4. The Extended Thiem’s Solution -Including the Impact of Heterogeneity think that an adaptive filter with high resolution, thus small filter volumes near the well and increasing filter volumes towards the far field, applies best to the flow equation. It still fulfills the physics of this strongly nonuniform flow pattern. As proposed by Attinger [2003] and Schneider and Attinger [2008] we use a Gaussian shaped filter function. Based on the nonuniform flow pattern for radial convergent flow we consider the filter length scale `to be proportional to the radius r`as described in Schneider and Attinger [2008], covering the idea of nearly no smoothing near the well and large smoothing in the far field. For isotropic media we choose a filter function of the form 2 fiso 1 x (x)= 3 exp -6, r 3/232 r2r with x 26R3, r`the radial distance from the well, treated as parameter and `a constant of proportionality, determining the width of the Gaussian filter. Note that fiso(x) is normalized. r In case of anisotropy we modify the filter in that way that in horizontal direction the filter width is 1,2 /6r, whereas in vertical direction the filter is weighted with the anisotropy ratio 3 /6e`·6r. At the same time the filter volume stays constant and independent of e, - thus 1 ·62 ·63 /6r3. Hence we find 1 = 2 = e` 13 r`and 3 = e` 23 r`which results in a more general filter function for anisotropic media 22-22 faniso 1 x1 + x2 + ex3 (x1,x2,x3)= 3 exp -6.`(4.7) r -2/32r2 3/23re What is the physical meaning of this spatial averaging? Imagine the filter in the original Coarse Graining procedure in three dimensions to be a block of volume 3 around a point x, where all spatial heterogeneities are averaged within this block, leading to a filtered head h(x). Then the radial depending filter can be seen as a cube around x growing with distance from the well. Near the well, the filter is very small leaving nearly all heterogeneity of K(x) unchanged, whereas far away the local conductivities are replaced by an averaged value. In the same way heads and fluxes are stepwise filtered with distance to the well. The influence of the anisotropy is covered by an adapted filter in vertical direction, giving that the filtering is still proportional to the correlation length in all directions. For e6 = 1 the fictitious cube filter around a point x is deformed to a cuboid with the same volume, reflecting the fact that the vertical compensation is reduced due to the stratification. Following the line of derivation in Attinger [2003] and Schneider and Attinger [2008] we evaluate KCG(r) with the correlation function in Eq.(4.3) and the adapted filter function in Eq.(4.7). A Details on the derivation are presented in the appendix A.1. We result in the filtered conductivity KCG(x)= KCG(r)+ K˜ CG(x,r), where K˜ CG(x,r) is the fluctuating part and KCG(r) is r AA the scale dependent mean hydraulic conductivity given by . . KCG A (r)= Kefu exp @2 (e) 1+ 3 . 22 !- r p6 32 2 `2 . ,`(4.8) e 4.4. The Extended Thiem’s Solution with the anisotropy function (e) and Kefu given in Eq.(4.5) and Eq.(4.4). 4.3.3. Representative Conductivity for Well Flow The Coarse Graining conductivity KCG(x)= KCG(r)+ K˜ CG(x,r) still contains local fluctu r A ations. Since we are interested in a representative conductivity for well flow which depends only on the radius, we perform a vertical average over a sufficiently thick aquifer. This reflects the fact that heads measured in real pumping tests can be regarded as means over the vertical extension of the well. Since the average of the fluctuating part K˜ CG(x,r) is zero by definition, KCG(r) remains as representative conductivity value for well flow. It can be seen in A Eq.(4.8) that the asymptotic behavior of KCG(r) covers KA as near field and Kefu as far field A representative conductivity values. As discussed in section 4.2.2 different boundary conditions at the well result in different near field representative conductivities. We therefore generalize the result to both possibilities of Kwell. For the BCA with Kwell = KA we use the abbreviation A = 2 (e) and for BCH with Kwell = KH we write H = 2( (e) - 1). We result in a general expression for the representative well flow conductivity for both boundary conditions . KCG(r)= Kefu exp @(2 ,e) 1+ 3 . 22 !- r v 32 . 2 . .`(4.9) e `2 A universal definition of `is given by `= ln(Kwell/Kefu). It allows KCG(r) to act as a general interpolating function between the near and far field representative conductivity values Kwell and Kefu, respectively. 4.4. The Extended Thiem’s Solution 4.4.1. Derivation of hefw(r) The Coarse Graining conductivity KCG(r), as given in Eq.(4.9) reflects the impact of heterogeneity on well flow in porous media. It is the representative conductivity value for which the solution of the flow equation best fits the drawdown of pumping tests. More precisely inserting KCG(r) to the flow equation (4.2) delivers a head which reproduces the vertically averaged drawdown of a three dimensional pumping test at every radial distance r`from the well in dependence on the parameters KG, 2,`and e`of K(x). The final step in our approach is to derive the effective well flow head hefw(r) by solving the radial flow equation with KCG(r). We transform Eq. (4.2) to polar coordinates, evaluate the vertical component and result in hefw(r) = C1 exp(-) ln r`+ C1 sinh()U1(r) (4.10) R` +C1 (1 - cosh()) U2(r)+ h(R) ,` 4. The Extended Thiem’s Solution -Including the Impact of Heterogeneity Figure 31: Comparison of hefw(r) and Thiem’s solution for two conductivity fields (KG = 10-4 m/s, 2 =1, e =1, f =5m, and f =10m) and both boundary conditions, BCA and BCH. Reference head is h(R =32m)=0m. with C1 =- Qw and 2LKefu r)+11 1 U1(r)=lnu(+ u(R)+1- u(r)u(R) , r)1 1 1 1 U2(r)=lnu( , u(R)- 2u(r)2 +2u(R)2 - 4u(r)4 +4u(R)4 q p where u(r)serves as abbreviation for u(r)=1+(r)2/(3e`)2 and . =ln(Kwell/Kefu). L is the aquifer thickness and Qw the pumping rate. h(R) is a reference head measured at the arbitrary distance R. This can be e.g. the head at the well h(R =rw)or a measured head value h(R)=hR in the far field, within the radius of influence of the pumping test. Details on the mathematical derivation can be found in the appendix A.3. The solution hefw(r)is a general expression for both boundary conditions at the well with hefw(r) A hefw for Kwell =KA (BCA) and hefw(r)for Kwell =KH (BCH). hefw(r)and (r)differ in the H AH parameter . For BCA we fix A =0.8. For BCH we set H =1.6. The relation A =0.5* H results from the fact that the transition zone from KA to Kefu is half the size of the transition from KH to Kefu. It results from the asymmetric relation between Kefu, KA and KH as discussed in section 4.2.2, see Figure 30. Analyzing hefw(r)in Eq.(4.10), we see that the first term resultsin Thiem’s Formula with Kwell as homogeneous substitute value. The terms U1 and U2 operate as correction terms, gaining in influence with increasing distance from the well. For homogeneous media . becomes zero and Eq.(4.10) reduces to Thiem’s solution, Eq.(4.1). 4.4. The Extended Thiem’s Solution Figure 32: Derivation of hefw(r)and hefw(r)with respect to parameters 2 and R in dimen- AH sionless scale r/`. The left plot shows the derivatives for 2 =0.5, the right one for 2 =2. Used setting: KG =10-4 m/s, e =1and h(r/R =6.4)=0m. As shown in Figure 31, hefw(r)interpolates between the drawdowns of Thiem’s solution with Kwell and Kefu as homogeneous substitute values, where the transition is determined by the correlation length `. In Figure 31 also the different sizes of the transition zone for BCA and BCH can be seen. 4.4.2. Sensitivity Analysis Formula (4.10) for the effective well flow head hefw(r) allows for a detailed analysis of the influence of the parameters of K(x) on the drawdown. Depending on Kwell meaning the applied boundary condition (BCA or BCH), the effects can be quite different. However, for both boundary conditions we state from our findings that the geometric mean KG influences the entire curve. The variance 2 mainly impacts on the drawdown at the well. And the correlation length R determines the ’velocity’ of transition from near to far field behavior. The variance 2 determines the magnitude of the drawdown in the vicinity of the well. For BCA a higher variance 2 causes larger values for KA and hence flattens the depression cone. For BCH the opposite effect appears. The larger 2, the steeper becomes the depression cone. This effect can be seen in Figure 32, where we plotted the absolute values of the derivative of hefw( r)with respect to the variance 2 for both near well conductivity values Kwell =KA and Kwell =KH and two choices of 2. For hefw(r), we can see that the influence of the variance A near the well is very strong and smooths out in the far field, but is still present through Kefu. For hefw H (r) the effect of 2 on the drawdown reverses because of the different signs in the exponent of KH and Kefu. In particular, within the distance of the first correlation length there exists an area where a change in 2 does not influence the depression of hefw(r)for BCH. The effect of a change in correlation length R on the drawdown can be seen in Figure 31 where we plotted hefw(r)for two different correlation lengths. The larger `, the longer takes the transition 4. The Extended Thiem’s Solution -Including the Impact of Heterogeneity Figure 33: Contour plot of anisotropy ratio e: the black lines show the isolines for the head hefw(r)in dependence on the dimensionless distance r/f and the anisotropy A ratio e. Used setting: KG =10-4 m/s, 2 =1and h(r/f =6.4)=0m. from Kwell to Kefu. In Figure 32 we plotted the dimensionless version of the derivative of hefw(r) with respect to the correlation length f for two different values of the variance 2. Noticeable is the increasing influence of f with increasing variance 2 for both boundary conditions BCA and BCH. This is caused by a larger distance between Kwell and Kefu for larger variances. Furthermore, it can be seen that the influence of f on hefw(r)vanishes quickly with increasing distance to the well. This is caused by the fact that the drawdown reaches the far field behavior after approximately two correlation lengths, meaning that hefw(r> 2`)= hThiem(r> 2`) with Kefu as homogeneous substitute value in Eq. (4.1). This is in line with the findings of Neuman et al. [2004, 2007]. If anisotropy (e6 =1) is assumed an additional quantity impacts the hydraulic head drawdown. The anisotropy rate e influences the far field behavior by its impact on Kefu in Eq.(4.4). Additionally, it is present as a scaling factor to the horizontal correlation length (visible in the expression of u(r)), since the relation between Kefu and Kwell also manipulates the transition zone. The impact of a change in e on the drawdown is plotted in Figure 33 for BCA, where we see that the sensitivity of hefw(r)towards e is very low. The same can be observed for BCH. A This can be explained by the fact that a stronger stratification does not impact very much on the flow pattern of pumping tests, which is mainly determined by horizontal flow. 4.4.3. The Head Solution as an Inverse Estimation Tool The analytical solution hefw(r)provides a useful tool to analyze pumping test data. Depending on the available amount of data, hefw(r)enables the inverse estimation of statistical properties under the assumption of a log-normal, Gaussian correlated hydraulic conductivity field. 4.5. Interpreting Numerical Pumping Test Examining Eq.(4.10) more closely we see that KG and 2 are both incorporated in Kwell  and Kefu, where the relation Kwell = KG exp±12 2only holds if ergodicity is fulfilled at the well. Generalizing hefw(r) to non-ergodic conditions, we shift the input variables from KG and 2 to Kwell and Kefu by using . = ln(Kwell/Kefu). However, from the discussion in section 4.4.2 it can be seen that Kwell, Kefu (respectively 2 and KG under ergodic conditions) and f all have a unique influence on the drawdown, which serves as a good basis for estimating them through a regression. Furthermore, Figure 32 shows that the estimation of f becomes even more certain the higher the variance 2 is. In contrast, e cannot be treated as an independent parameter, because of its low influence on hefw(r), see Figure 33. We therefore support the statement of Firmani et al. [2006] that an estimation of the anisotropy ratio e through pumping test data is very error prone. What is not discussed until now is the influence of the boundary condition assumed at the well on the parameter estimation. The certainty in estimating Kwell and Kefu (respectively KG and 2) is independent on the choice of BCA or BCH. But to infer on f the crucial area is the transition zone from Kwell to Kefu. This zone is smaller for BCA than for BCH, independent of 2, because of the relation of Kefu to KA and KH, see Figure 30. In case of anisotropy the transition zone for BCA becomes even smaller and vanishes for stratified media since Kefu(e = 0) = KA. This makes a parameter estimation with BCA more difficult and less reliable. For BCH the opposite effect occurs. A convergence of Kefu to KA enlarges the transition zone, thus improves the ability of estimating `. 4.5. Interpreting Numerical Pumping Test In this section we analyze three dimensional numerical pumping tests in highly heterogeneous porous media to confirm the validity of hefw(r) to describe well flow effectively. We examine the drawdown results and develop a method to infer on the statistical parameters of the conductivity distribution K(x). With our simulations we characterize the transition zone of near well to far field behavior being the main area of influence of heterogeneity. Furthermore, we discuss the question if a single pumping test realization is sufficient to infer on all parameters of K(x) and present several ways to cope with a lack of ergodicity. 4.5.1. Simulation Setup To examine the question of ergodicity and the influence of domain size on a numerical pumping test drawdown, we perform simulations on several meshes of different horizontal and vertical extension. Starting from the findings of Firmani et al. [2006], we generate a large domain G3 with horizontal mesh size of R3 = 40`, a medium sized domain G2 with R2 = 25.6f and a small one, G1 with R1 =6.4`, all of them having a vertical extension of 64`v. Adapting to the radial flow system, we establish a mesh refinement in the range of -f to `. It provides a high resolution of the head drawdown near the well. Additionally, the well is not included as a point source but as a hollow cylinder with radius rw =0.01m. All meshes refer to a correlation 4. The Extended Thiem’s Solution -Including the Impact of Heterogeneity Ensemble KG [10-4 m/s] 2 i [m] e E1 1.0 1.0 5.0 1 E2 1.0 2.0 5.0 1 E3 1.0 1.0 10.0 1 E4 1.0 1.0 10.0 0.5 Table 17: Input parameters of generated ensembles. length of i = 5 m and a resolution of 5 cells per correlation length, resulting in a cell size of 1m, being in the range of an idealized small scale pumping test. Note that for the ensembles with a correlation length of i = 10 m these ratios change while keeping the meshes unchanged. The boundary conditions of the simulations are the following: The upper and lower horizontal planes delimiting the domain are set impervious which reflects a confined aquifer. At the radial distance Ri according to Gi a constant head h(Ri) = 0 is applied, giving a circular outer horizontal boundary condition. At the well we use a constant total pumping rate of Qw = -10-3 m3/s with two different boundary conditions: (i) constant flux (BCH), where we assign the same pumping rate Qi at every grid cell, resulting in Kwell = KH and (ii) proportional flux (BCA), where the assigned rate is proportional to the local conductivity of the grid cell Qi . Ki, giving Kwell = KA, as stated by Indelman et al. [1996] and Firmani et al. [2006]. All simulations are performed using the finite element software OpenGeoSys developed by Kolditz et al. [2012a]. The code was tested against a steady state pumping test with homogeneous conductivity and the results in two and three dimension are in very good agreement with the analytical solution of Thiem (4.1). To generate heterogeneous, log-normal distributed, Gaussian correlated conductivity fields we make use of the statistical field generator randomfield provided by Cirpka [2010]. Hydraulic conductivity fields are created by a given deterministic power spectra, as described in Dykaar and Kitanidis [1992a]. We generate several realizations of hydraulic conductivity fields in three dimensions of the same statistical parameters, i.e. geometric mean KG, variance 2, correlation length i and anisotropy ratio e all forming one ensemble. To investigate the influence of parameters on the drawdown numerically, we generate several ensembles with varying parameter setups, listed in Table 17. In this study every ensemble consists of 20 realizations. RL We focus on the vertical average of the simulated head hh(x1,x2)) =1/L0 h(x1,x2,x3)dx3. We take the mean of hh(x1,x2)) on the x1 and x2 axes as representative drawdown hh(r)) for a realization, only depending on the distance to the well r. In a first step, we compare the drawdown hh(r)) with hefw(r) and Thiem’s solution hThiem(r), using Kefu and Kwell as homogeneous substitutes. In a second step, we use hefw(r) to infer on the statistical parameters of K(x) by applying a nonlinear regression to find estimates Kˆ G,ˆ2 ,`ˆwhich minimize the mean square error between the numerical drawdown data hh(r)) and hefw(r). In this estimation procedure we include all available points ri in the range of 5i (corresponds to 25m) of the drawdown hh(r)i. The question of the applicability of hefw(r) on limited head data is of quite a complex nature. Answering it from the perspective of numerical pumping tests 4.5. Interpreting Numerical Pumping Test Figure 34: Comparison of 20 simulated drawdowns (BCH on G1) from ensemble E2 (KG = 10-4 m/s, 2 =2, g =5m, e =1) in log scale. with full data availability results in a complicated selection criteria. The necessary detailed statistical analysis is given in section 5. Though in reality a full range of head data will not be available, we first focus on examining the validity of hefw(r) more theoretically in order to describe three dimensional well flow effectively. 4.5.2. Influence of Domain size We test different horizontally extended domains to investigate how the position of the outer boundary condition influences the simulated head. Comparing the numerical drawdown results for identical heterogeneity fields on all three domains G1, G2 and G3 shows negligible differ ences in the qualitative behavior, valid for all realizations and all ensembles. A quantitative hhG1(r)i-hhG3 (r)) between the hhG1(r)i comparison proves that the maximum relative difference maxr drawdowns on G1 and G3 is less than 0.5%. Furthermore, the inverse estimation results are nearly identical (not shown here). The numerical simulations -independent of the used domain -confirmed the findings of the sensitivity analysis in 4.4.2 that the mean drawdown reaches the far field representative behavior after less than two correlation lengths. It shows that reducing the horizontal extension does not impact on the simulated drawdown, as long as the outer boundary applies at a distance of more than tree to four correlation length. In the following we will thus reduce the discussion on numerical results to domain setup G1, since we believe that it contains all information necessary to determine the impact of heterogeneity on the drawdown. In the next step, we test the influence of the vertical domain extent on the drawdown results. We observe significant differences in the mean drawdown when reducing the vertical domain size, being in the line with the findings of Firmani et al. [2006]. They state that a vertical extension 4. The Extended Thiem’s Solution -Including the Impact of Heterogeneity Figure 35: Plot of simulated drawdowns (BCA and BCH on G1) of a single realization of ensemble E1 (KG = 10-4 m/s, 2 =1, f =5m, e =1) versus hefw(r) with local and theoretical value for near well conductivity: hKAi=1.422 ·10-4 m/s, KA = 1.649 ·10-4 m/s, hKHi=0.523 ·10-4 m/s, KH =0.607 ·10-4 m/s. of L 60`z ensures ergodicity. To check whether this also holds true for our simulation setup with L being even 64`z we compare the resulting drawdowns hh(r)ifor one ensemble and find large differences between the 20 realizations. We find that behavior in particular near the well, as shown in Figure 34 for ensemble E2, which is in contradiction to the findings of Firmani et al. [2006]. The reason for the spreading within one ensemble becomes evident when relating the drawdown hh(r)iof every realization to its local hydraulic conductivity value at the well hKwelli. P 1 m According to both boundary conditions, we calculate the arithmetic mean hKAi= mi=1 Ki -1 P 1 m 1 and the harmonic mean hKHi=i=1 of the conductivity values Ki of the m = 320 mKi cells along the well. As predicted by theory, the local mean hKwellidetermines the depression in the vicinity of the well. But in contrast to Firmani et al. [2006], who found differences of  less than 1.4% between hKwelliand the theoretical value Kwell = KG exp± 12 2 , we found differences up to 20%.  The discrepancies between local and theoretical expected means, hKwelli6= KG exp± 12 2 , we trace back on a too small sample size. A mean over 320 K-values, moreover spatially correlated, is not representative for a full K-field of over one million values, ensuring the convergence to the theoretical expected mean. We therefore state that a single three dimensional pumping test, with a randomly chosen well position, does not fulfill ergodic conditions, even for very large vertical domain extensions. 4.5. Interpreting Numerical Pumping Test Figure 36: Statistics on results of ensemble E1 (n = 20): The upper boxplots compare the local means at the well hKwell) with the inverse estimated values Kˆ well. The lower boxplots compare the estimation results of Kˆ efu and `ˆfor both boundary conditions (BCA and BCH). The appropriate theoretical value is marked on the vertical axis. 4.5.3. Analysis of a Single Pumping Test To overcome the lack of ergodicity, we analyze single pumping tests with hefw(r), accounting for local variations at the well. We use the universal definition of . = ln(hKwelli/Kefu) as introduced and discussed in section 4.3.3 and 4.4.3 with the local mean hKwell) instead of the theoretical value Kwell and find very high accordance between hh(r)) and hefw(r). Figure 35 shows the impact of hKwell) on the simulated drawdown of a single realization for ensemble E1 and the large differences between hefw(r) when using hKwell) compared to Kwell. It can be seen quite well that hefw(r) with hKwell) matches the simulated drawdown for both boundary conditions. This is clearly not the case for hefw(r) with Kwell, because of the significant differences of 13.8% between Kwell and hKwell) for both boundary conditions. It should be mentioned that a modification of hefw(r) goes along with a shift of input parameters from KG and 2 to the far and near field representative values hKwell) and Kefu. In terms of analyzing single realizations, it can be interpreted as a decoupling of the full field variance 2 incorporated in Kefu and a local variance at the well h2) incorporated in hKwelli. Both values, local and full field variance, might differ significantly from each other. We performed a parameter estimation by using Eq. (4.10) with . = ln hKwell) on every realiza- Kefu tion of ensemble E1. A statistical analysis of the results is shown in the box plots of Figure 36. For both boundary conditions, we find a very high accordance between the local mean values at the well hKwell) (hKA) and hKHi, respectively) and the estimated values Kˆ well (Kˆ A and Kˆ H), plotted in the upper box plots. The estimation results for the far field conductivity Kˆ efu are 4. The Extended Thiem’s Solution -Including the Impact of Heterogeneity ˆKG [10-4 m/s] ˆ2 ˆKefu [10-4 m/s] ˆKA [10-4 m/s] `ˆ [m] E1 1.0 1.0 1.181 1.649 5.0 0.972 1.064 1.16 1.654 4.99 (±0.014) (±0.028) (±0.013) (±0.008) (±0.42) E2 1.0 2.0 1.396 2.718 5.0 0.946 2.083 1.338 2.68 5.05 (±0.023) (±0.05) (±0.023) (±0.028) (±0.41) E3 1.0 1.0 1.181 1.649 10.0 0.954 1.079 1.142 1.636 10.12 (±0.005) (±0.011) (±0.004) (±0.005) (±0.55) Table 18: Parameter estimation results of simulated mean head hh(r)) for three ensembles with BCA. Expected ensemble parameters in italic, the inverse estimates in bold and the 95% confidence intervalls in brackets. shown in the lower left box plot. We can see that BCA underestimates and BCH overestimates the theoretical expected value Kefu, but only by small deviations. The estimated values for the correlation length `ˆ match the theoretical value f very well, visible in the lower right box plot. However, the variability is quite large, especially for BCA. This corresponds to the large spread of hKA) from KA, because the estimation of the correlation length depends strongly on the transition zone and is therefore triggered by the local discrepancies at the well. We conclude from the analysis of a single realization that hefw(r) allows a good estimation of the statistical parameters Kˆ well, Kˆ efu and `ˆ. Since Kˆ well is very much influenced by the local distribution of the conductivity near the well, it does not necessarily correspond to the  theoretical value Kwell = KG exp±21 2and we thus recommend to be careful when tracing back KG and 2 from Kˆ well and Kˆ efu for a single pumping test. 4.5.4. Analysis of an Ensemble of Pumping Tests Since a vertical extension of L = 64`z does not ensure ergodicity sufficiently, we state that even in three dimensions it is necessary to investigate an ensemble of pumping tests to infer on the statistics of K(x). A further extension in vertical direction would also be a possibility to ameliorate the results in theory. However, even assuming high anisotropy rates these conditions will hardly be fulfilled in reality. On the other hand, a number of pumping tests within an observation area will give much better insights into the heterogeneous structure of the subsurface. Within our setting we tested several ensemble sizes up to n = 20 realizations and estimate P n the statistical parameters for the ensemble mean hh(r)i=i=1hhi(r)i. We found a quick n convergence of Kˆ efu and `ˆto the theoretical assigned values: less than 10 realizations are sufficient to be in a range of 95% accuracy. For the near well conductivity about 15-20 realizations are necessary to ensure an acceptable good agreement between Kˆ well and Kwell. Table 18 and Table 19 show results for the estimated parameters Kˆ well, Kˆ efu and `ˆfrom hh(r)in=20 for three different ensembles and both boundary conditions. The values for Kˆ G and ˆ2 can either be 4.6. Summary and Conclusions ˆKG [10-4 m/s] ˆ2 ˆKefu [10-4 m/s] ˆKH [10-4 m/s] `ˆ [m] E1 1.0 1.0 1.181 0.607 5.0 1.008 0.972 1.186 0.62 5.01 (±0.003) (±0.005) (±0.004) (±0.001) (±0.09) E2 1.0 2.0 1.396 0.368 5.0 0.986 2.025 1.383 0.358 5.07 (±0.007) (±0.013) (±0.012) (±0.001) (±0.08) E3 1.0 1.0 1.181 0.607 10.0 1.031 0.94 1.206 0.644 10.26 (±0.006) (±0.011) (±0.009) (±0.001) (±0.27) Table 19: Parameter estimation results of simulated mean head hh(r)) for three ensembles with BCH. Expected ensemble parameters in italic, the inverse estimates in bold and the 95% confidence intervalls in brackets. evaluated from Kˆ well and Kˆ efu or directly estimated by using A =  (e) and H = ( (e)-1); the estimation results are also listed in Table 18 and Table 19. Again we state that the good estimation results for Kˆ well and ˆ2 are mainly caused by the fact that for a sufficiently large number of realizations, the local mean at the well hKwell) converges to the theoretical assigned value Kwell. A precise number of pumping tests needed to ensure ergodicity can be traced back on the question: What sample size N is necessary to guarantee the convergence of the mean of a sample of spatially correlated log-normal distributed  1 PN 1 values hKiN = i=1 Ki to the theoretical expected mean KA = KG exp2 2. Answering N this question is out of the scope of this study. Most likely, there does not exist a single number being appropriate for all possible statistical and geometrical settings. An item which is not discussed until now is the interpretation of drawdowns in anisotropic media. As discussed in section 4.4.2 we do not incorporate e into our estimation procedure due to the low sensitivity of hefw(r) towards e. However, applying hefw(r) on drawdown data in anisotropic media is possible and useful, as shown in Figure 37. Assuming a reasonable value for e leads to very good accordance of hh(r)in=20 and hefw(r) and further allows the estimation of Kˆ efu, Kˆ well and `ˆ. For application on real pumping tests data we would recommend to fix e to a reasonable value or to carry out several estimations with various ratios for e like 0.01, 0.1, 0.5 and 1.0 and interpret the results with respect to the accordance of Kˆ efu, Kˆ well and `ˆ to the drawdown. We conclude by stating that our numerical results show that hefw(r) is a promising tool to characterize aquifer properties like mean conductivity, variance, and spatial correlation at a very local scale by interpreting the near well behavior of steady state pumping tests. 4.6. Summary and Conclusions In this study we introduced a representative description of the hydraulic head drawdown for a steady state pumping test with fully penetrating well for highly heterogeneous media. By making use of the upscaling technique Coarse Graining, we derived a radial depending con 4. The Extended Thiem’s Solution -Including the Impact of Heterogeneity Figure 37: Plot of ensemble drawdown hh(r)in=20 for ensemble E4 (KG = 10-4 m/s, 2 = 1, f = 10 m, e =0.5) versus hefw(r) in anisotropic media for BCA and BCH. The inset shows the near well behavior in log scale. ductivity KCG(r). It interpolates between the known near and far field representative conductivities for well flow. Therefore, we deduced the effective well flow head solution hefw(r) which reproduces the mean drawdown of a pumping test adequately. We understand hefw(r) as an extension of Thiem’s Formula incorporating the effects of the statistical parameters of the underlying log-normal distributed conductivity field K(x) on the flow pattern. The analytical character of hefw(r) allowed us to perform a sensitivity analysis for the parameters of K(x) on the drawdown. We found that the variance 2 has the strongest impact on the hydraulic head directly at the well. The horizontal correlation length f determines the transition from near to far field behavior. In particular the impact of f increases with increasing variance 2, which makes a prediction of f easier for highly heterogeneous media. The anisotropy ratio e has only little influence on the drawdown, giving that hefw(r) shows very low sensitivity towards changes in e. To validate the applicability of hefw(r) we performed steady state numerical pumping tests in three dimensional highly heterogeneous anisotropic media, with variances up to 2 =2. Our investigations confirmed the findings of Indelman et al. [1996] that the far field behavior  is covered by Kefu = KG exp2( 12 - (e)). We also found the near well representative conductivity Kwell to be the arithmetic KA or harmonic mean KH, depending on the assigned Dirichlet or Neumann boundary condition. However, the means at the well have to be considered very locally. Investigations on the local distribution of K(x) showed that the arithmetic and harmonic mean of the conductivity values directly along the well hKA) and hKH) are not representative for the theoretically expected  means of the full field KA = KG exp 12 2 and KH = KG exp - 12 2 . We found discrepancies up to 20% between KA and hKA) as well as between KH and hKH) for all tested variances. We therefore conclude that a single pumping test realization does not fulfill ergodic conditions in the vicinity of the well even for a large vertical extension of more than 60 correlation lengths. 4.6. Summary and Conclusions This is in contrast to a previous work published by Firmani et al. [2006]. In order to make predictions for the overall statistics, we analyzed ensembles of pumping tests and showed that hefw(r) does not only reproduce the ensemble drawdown but enables the estimation of the statistical parameters with very high accuracy. In isotropic media the estimated results Kˆ G, ˆ2 and `ˆ differ less than 6% from the expected ones for all ensembles with very high confidence. Solely the anisotropy ratio e is difficult to infer. We agree with Firmani et al. [2006] that the estimation of e from pumping test data is very error prone. Nonetheless, hefw(r) also allows the interpretation of drawdowns in anisotropic media (e< 1) by assuming a reasonable ratio e and then estimating the parameters Kˆ G, ˆ2 and `ˆ. However, being limited to ensemble averages of multiple pumping tests is clearly a limitation for interpreting real drawdown data. To overcome the lack of ergodicity at the well in a single realization we adapted our proposed formula hefw(r) on local statistics of the conductivity, incorporating hKA) and hKHi, which gives a much better reproduction of the simulated depression cone than with theoretically values KA and KH. This modification allows to estimate Kˆ A and Kˆ H respectively, Kˆ efu and `ˆ for single drawdown data. If several pumping tests in one area are available each can be interpreted with hefw(r) and afterwards a statistical analysis can be applied to infer on ˆ2 and Kˆ G. Thus hefw(r) can serve as a helpful tool to interpret real drawdown data for an arbitrary number of steady state pumping tests. Exploiting our results with respect to predictions on a real pumping test sampling design we suppose that the quality of the parameter estimation mainly depend on the position of the observation wells. A good estimation of the variance 2 requires measurements directly at the well. To infer on the correlation length f the vicinity of the well, meaning the area within two correlation lengths has to be investigated. Measurements far from the well allow to infer on Kefu. The larger the number of head data in the corresponding area of influence of a parameter, the more reliable are its estimation result. Thus we can use hefw(r) not only to infer on the statistics but it also allows to judge the usefulness of measurements with respect to the estimation of the parameters for the underlying hydraulic conductivity field. 4. The Extended Thiem’s Solution -Including the Impact of Heterogeneity 5. Estimating Parameters of Aquifer Heterogeneity Using Pumping Tests -a Paradigm for Field Applications 5.1. Introduction Groundwater pumping tests are applied as common and well-established tool to estimate hydraulic properties of porous media. Estimating parameters such as the saturated hydraulic conductivity of an aquifer is based on measurements of the hydraulic head across the area of drawdown. Classical approaches to estimate hydraulic conductivity from pumping tests are based on Thiem’s formula [Thiem, 1906] or Theis’ solution [Theis, 1935] for steady state or transient flow conditions, respectively; the latter having been simplified by Cooper and Jacob [1946]. Detailed information on the application of those approaches under numerous boundary conditions was provided by Kruseman and de Ridder [2000]. With regard to natural aquifers, a critical shortcoming of these methods is the simplified assumption of aquifers being homogeneous. The vast majority of natural aquifers are characterized by discontinuities evolved from geomorphologic processes which affect the drawdown curve derived from pumping tests [Dagan and Neuman, 1997]. Estimating the heterogeneous structure of aquifers is crucial for characterizing the groundwater flow toward the pumping well and for predicting solute transport, particularly in regard to contaminant migration and in-situ remediation. For example, assuming a log-normal hydraulic conductivity model, the horizontal dispersivity of an aquifer is proportional to the product of the variance 2 and the horizontal correlation length f [Gelhar, 1993]. Taking heterogeneity into account, the hydraulic conductivity is assumed to be represented by a spatial random function K(x) with log-normal distributed values K(x) . LN (µ, 2), where µ and 2 denote the mean and variance of the normal distribution ln K(x), respectively. However, using a single representative mean value of K(x) to describe the water flow toward the pumping well is inadequate for steady state conditions [Matheron, 1967], because distinct representative values emerge for the flow near (Kwell) and far from the pumping well (Kefu). The mathematical expressions of Kwell and Kefu depend on µ and 2 and are related to the dimensionality of the water flow, see Sánchez-Vila et al. [2006] for details. Most studies concerning well flow in heterogeneous porous media are limited to a two dimensional flow model and therefore are representative for large scale pumping tests, e.g. Desbarats [1992]; Sánchez-Vila et al. [1999]; Copty and Findikakis [2004]; Neuman et al. [2004]; Leven and Dietrich [2006]; Neuman et al. [2007]; Dagan and Lessoff [2007]; Schneider and Attinger [2008]. However, for small scale pumping tests vertical flow in the vicinity of the pumping well is a critical component that influences the drawdown. Those considerations indicate the need to reliably estimate statistical parameters of K(x) using data from pumping tests and applying simple analytical solutions for three dimensional water flow [Indelman and Abramovich, 1994; Indelman et al., 1996; Indelman, 2001; Guadagnini et al., 2003; Firmani et al., 2006; Zech et al., 2012], thereby reducing time and cost-load otherwise involved with methods such as laboratory trials. 5. Estimating Parameters of Aquifer Heterogeneity Using Pumping Tests By deriving the vertical mean drawdown of a steady state pumping test in three dimensional heterogeneous porous media, we presented in the previous section 4.4 an analytical tool to analyze statistical parameters of K(x). The developed formula – referring to the effective well flow head hefw(r) – depends on the radial distance from the pumping well r and the statistical properties KG, 2, c and e of K(x). In contrast to Indelman and Abramovich [1994] and subsequent works, this method is based on an upscaling procedure called Coarse Graining [Attinger, 2003] rather than perturbation theory. Numerical simulations of pumping tests showed that hefw(r) reliably reproduces the radial depression cone of a steady state pumping test for a range of statistical parameter values. In contrast to simulated pumping tests, the underlying K(x) distribution of on site aquifers is unknown, which hampers the qualitative assessment of parameter estimates when using empirical field-site data. Moreover, sample size of head measurements is limited, which involves uncertainty in parameter estimates, e.g. if pumping tests are by chance conducted in areas of extraordinary low or high permeability. Thus the question remains whether hefw(r) is capable to estimate the statistical parameters of an natural aquifer using field-site pumping test data. The aim of this study is to close the gap between theoretical and field application of hefw(r) on pumping test data. We examined the capability and predictive power of hefw(r) to provide estimates of statistical parameters of K(x) from pumping test data such as under field-site conditions. As a first step (section 5.3), we analyzed simulated pumping tests by reducing the sample size of head measurements and evaluating the quality of parameter estimation. As a second step (section 5.4), we applied hefw(r) on empirical groundwater pumping test data from the field site Horkheimer Insel, Germany [Schad and Teutsch, 1994] and compared the hefw(r) estimates with the findings from laboratory investigations [Schad, 1997]. Based on these findings, we provide consideration regarding the conceptual design of groundwater pumping tests and the predictive power of established pumping tests sites. 5.2. Theoretical Framework This study is based on the effective well flow head hefw(r) introduced in section 4 and published in Zech et al. [2012]. The relation between the effective well flow head and the hydraulic head for pumping tests conducted in heterogeneous porous media is depicted in Figure 38. Moreover, the fundamentally distinct method of predicting the hydraulic head of homogeneous aquifers using Thiem’s solution is illustrated. Starting point is a spatially distributed heterogeneous conductivity field K(x,KG,2, `, e). Conducting pumping tests on such an aquifer results in a spatially distributed hydraulic head field h(r, , KG,2, `, e), also depending on the statistics of K(x). Tracing back the parameters of K(x) from h(r, , KG,2, `, e) is the fundamental goal of estimating aquifer parameters inversely. The approach of Zech et al. [2012] is based on adapted spatial averaging of K(x,KG,2, `, e) according to the conditions of pumping tests. Using the upscaling procedure Coarse Graining [Attinger, 2003] a representative conductivity KCG(r, KG,2, `, e) is derived, which is not fully homogenized, but still depends on the statistical parameter of K(x). 5.2. Theoretical Framework Figure 38: Relation between the deduction of hefw(r)and processes occurring during field site pumping tests. Coarse Graining can be best explained as a spatial filtering procedure that averages over volumes of variable filter size. Applied to pumping tests, Coarse Graining allows to account for the character of radial convergent flow. Near the pumping well, where the impact of heterogeneity is large, the filter size is adjusted such as that small volumes are captured, thereby leaving the heterogeneity unchanged. Whereas far from the pumping well, the filter size is larger resulting in averaged values. Performing a pumping test on KCG(r)results in the effective well flow head hefw(r, KG,2, `, e) (Figure 38). Numerical simulation in section 4.5 showed that hefw reproduces the vertical averaged hydraulic head h(r, , KG,2, `, e)of a pumping test in randomly heterogeneous medium, due to its ability to capture effects of heterogeneity in contrast to Thiem’s solution hT(r, KG). 5.2.1. The Effective Well Flow Head The effective well flow head hefw(r) can be considered an extension of Thiem’s solution to heterogeneous media. It reproduces the vertical mean drawdown of a steady state pumping test in relation to the radial distance r and the statistical parameters of K(x): hefw( r)=C1 exp(-)ln r +C1 sinh()U1(r)+C1 (1- cosh())U2(r)+h(R) , (5.1) R r Qw , . =lnKwell r with the abbreviations C1 =- , u(r)=1+pand the terms 3 2LKefu Kefu e`2 r)+11 1 U1(r)=ln u(+ u(R)+1- u(r)u(R) r)1 1 1 1 U2(r)=ln u( , u(R)- 2u(r)2 +2u(R)2 - 4u(r)4 +4u(R)4 5. Estimating Parameters of Aquifer Heterogeneity Using Pumping Tests where Qw is the pumping rate, L is the aquifer thickness, . is a constant of proportionality, f is the horizontal correlation length and e . [0, 1] is the anisotropy ratio assuming a Gaussian shaped spatial correlation structure in horizontally isotropic and vertically anisotropic media. The reference head h(R) is measured at the radial distance R, which can either be at the pumping well or in the far field. The hefw(r)-solution interpolates between the representative hydraulic conductivities at the pumping well and in the far field – Kwell and Kefu, respectively. They differ remarkably for steady state flow conditions, thereby affecting the flow characteristics of each radial section of the depression cone. For log-normal distributed conductivity the asymptotic drawdown behavior in the far field of a steady state pumping test in three dimensions is characterized by the effective hydraulic conductivity for uniform flow [Dagan, 1989; Indelman and Abramovich, 1994; Sánchez-Vila et al., 2006]:  1 Kefu = KG exp 22 - (e), (5.2) where KG denotes the geometric mean, 2 the variance and (e) is the anisotropy function p e v 1 (e)= 2 arctan1/e2 - 1- e. 2(1-e2)1-e Estimates of Kwell depend on the assigned boundary condition, [Indelman et al., 1996; Indelman and Dagan, 2004]: (i) the Neumann boundary condition (hereafter referred to as BCH) is based on constant flux and corresponds to the harmonic mean Kwell = KH, whereas (ii) the Dirichlet boundary condition (BCA) is based on a constant head and corresponds to the arithmetic mean Kwell = KA. If ergodicity applies at the pumping well the values for Kwell can be  calculated as KH = KG exp - 12 2 for BCH, and KA = KG exp 12 2 for BCA. Both bound ary conditions are addressed in Eq.(5.1) through Kwell in . = ln Kwell and the parameter , Kefu which is BCH =1.6 and BCA =0.8, see section 4.4.1. Notably, hefw(r) is a deterministic head solution governed by the statistical characteristics of K(x) rather than a spatial random function. The solution presented in Eq.(5.1) is based on the four parameters Kefu and Kwell, f and e. If ergodic conditions apply 2 and KG can be calculated via 2 = 2 (ln Kwell - ln Kefu) 2 (e) - 1 ± 1 , (5.3)  with + for BCA and - for BCH and KG = Kwell exp = 12 2-2 = Kefu exp 12 - (e) . However, for single pumping tests ergodicity cannot be presumed in the vicinity of the pumping well and Kefu and Kwell are required for Eq.(5.1) to be independent from any ergodicity assumption. This modification is particularly useful for single pumping tests where Kwell is not related to the theoretical values of KH or KA (section 4.5.2). Both parameterizations, KG and 2 for ergodic conditions, and Kefu and Kwell for non-ergodic conditions, are able to reproduce ensemble head means of simulated pumping tests. Since the aim of this study is to examine the capability and predictive power of hefw(r) to estimate statistical parameters using single pumping test data from field campaigns, we focus on non-ergodic conditions hereafter. 5.3. Pumping Tests under Limited Data Availability -a Paradigm for Field Applications 113 5.2.2. Inverse Parameter Estimation While the analytical character of hefw(r) allows for examining the influence of statistical parameters of K(x) on the depression cone, it also can be utilized to estimate those parameters inversely from pumping tests. Based on this premise nonlinear regression is used to estimate best values of Kˆ efu,Kˆ well and `ˆto minimize the mean square error of the difference between hefw(r) and the measured drawdown hh(r)i. However, to determine the accuracy and certainty of those estimates it is essential to ascertain the sensitivity of hefw(r) to each of the parameters. A sensitivity analysis, performed in section 4.4.2, revealed that Kwell and Kefu influence hefw(r) in the vicinity of the pumping well and the far field, respectively, whereas the correlation length f controls the transition between Kwell-driven and Kefu-driven drawdown. The larger f the larger is the transition zone around the pumping well. Moreover, hefw(r) is highly sensitive to changes in values of Kefu, Kwell and `, whereas the sensitivity of hefw(r) towards changes in the anisotropy ratio e is very low. Hence, estimating e from pumping tests is hardly possible [Firmani et al., 2006; Zech et al., 2012] and reliable estimates of Kˆ efu, Kˆ well and `ˆ can only be achieved if e is considered constant. Comparing the confidence intervals, denoted by CI, of the three parameters reveals that CI`ˆ is much greater than CI ˆand CI ˆ. This is caused by the distinct areas that effect Kˆ efu Kefu Kwell and Kˆ well – far field and vicinity of pumping well, respectively – whereas the accuracy of `ˆalso depends on the estimates of Kˆ efu and Kˆ well. Alternatively, to eliminate this intrinsic uncertainty and thereby decreasing the confidence interval, the correlation length `˜ rather than f can be estimated using Kˆ efu and Kˆ well as constant parameters. 5.3. Pumping Tests under Limited Data Availability -a Paradigm for Field Applications We examined the capability of hefw(r) to estimate statistical parameters of a heterogeneity under limited data availability as found under conditions of on-site pumping tests. In order to characterize the impact of the number and location of head measurements on the estimation quality, represented by accuracy and uncertainty, we performed simulated numerical experiments. The simulations were based on head data of pumping tests conducted in three dimensional aquifers with randomly generated hydraulic conductivity fields. Sample size of head measurements was systematically reduced with regard to their horizontal distribution and their radial distance relative to the pumping well. 5.3.1. Pumping Test Model The finite element software OpenGeoSys [Kolditz et al., 2012a] was used to simulate pumping tests on randomly generated realizations of three dimensional hydraulic conductivity fields. The -1 statistical parameters of the aquifer were set as KG = 10-4ms, 2 = 1, f =8m, and e = 1. 5. Estimating Parameters of Aquifer Heterogeneity Using Pumping Tests r i in[m] #r i characteristics S1 r w, 1,..,32 33 full range, equidistant, 1m interval S2 r w, 2,..,32 17 full range, equidistant, 2m interval S3 r w, 4,..,32 9 full range, equidistant, 4m interval S4 r w, 8,..,32 5 full range, equidistant, 8m interval S5 r w, 1,..,8 9 vicinity of the well S6 9,..,32 24 far field S7 r w, 2, 4, 6, 8, 16, 24, 32 8 full range, focus on vicinity of the well S8 r w, 8, 16, 20, 24, 28, 32 7 full range, focus on far field Table 20: Scenarios S1-S8 of limited radial head data r i. The total number of measurements is denoted as #r i. Head measurements were sampled from all four axial directions. The numerical grid was similar to the one used previously in section 4.5, with small adaption to conditions of short term pumping tests as conducted on-site at the Horkheimer Insel by Schad [1997], (section 5.4). The horizontal grid size was 8H (64m) with a uniform grid cell size of 18 H (1m) except for cells in the vicinity of the pumping well, which were required to be smaller for the purpose of integrating steep head gradients. The vertical grid size was 15eH with a uniform grid cell size of 18 eH . Impermeable horizontal layers formed the base and top of the aquifer. The outer boundary condition was set as fixed head h (R =4H ) = 0. At the pumping well both BCH and BCA 3-1 were applied, respectively (section 5.2). For the first the pumping rate of Q w = -10-3ms was equally distributed over all vertical grid cells at the pumping well, whereas for the latter the pumping rate was proportional to the vertically distributed hydraulic conductivity values at the pumping well Q i . K i, which corresponds to the constant head boundary condition [Firmani et al., 2006; Zech et al., 2012]. Simulated heads were measured at the four axial directions . 1 =0,. 2 =1/ 2, . 3 = , . 4 =3/ 2p with a resolution of 1m. We used vertically averaged heads to reflect on head measurements as derived under field-site conditions, eventually only depending on the horizontal position de RL scribed by the radial distance r i and the angular position . j: hh (r i,. j)) =1/L0 h (r i,. j,z )dz . 5.3.2. Methods of Sample Size Reduction We reduced the sample size of simulated head measurements stepwise according to (i) their angular position, (ii) radial distance relative to the pumping well, and (iii) a combination of both. For each step we estimated the statistical parameters Kˆ efu, Kˆ well and `ˆaccording to the procedure described in section 5.2.2. The average hydraulic head hh (r )) derived as the mean of the four axial directions . 1 = 0, . 2 =1/ 2p , . 3 = p , . 4 =3/ 2p and measured at 1m intervals serves as reference scenario S1 (Table 20) with full data availability. Firstly, we controlled the sample size of head measurements with respect to their angular arrangement, while remaining the full sample size with regard to their radial distance from the pumping well. The statistical parameters were estimated for each of the four axial directions, 5.3. Pumping Tests under Limited Data Availability -a Paradigm for Field Applications 115 ˆKwell [10-4 m/s] ˆKefu [10-4 m/s] `ˆ [m] reference parameter 0.648 1.181 8.0 hh(r)) 0.665 (±0.003) 1.165 (±0.012) 8.86 (±0.54) hh(r, 1)) hh(r, 2)) hh(r, 3)) hh(r, 4)) 0.657 (±0.006) 0.663 (±0.006) 0.677 (±0.009) 0.664 (±0.005) 1.176 (±0.020) 1.207 (±0.027) 1.199 (±0.054) 1.106 (±0.016) 7.58 (±0.83) 8.95 (±1.07) 12.52 (±2.42) 7.60 (±0.82) h˜h(r, C1)) h˜h(r, C2)) h˜h(r, C3)) h˜h(r, C4)) h˜h(r, C5)) 0.660 (±0.011) 0.659 (±0.010) 0.661 (±0.012) 0.672 (±0.011) 0.660 (±0.010) 1.155 (±0.039) 1.152 (±0.034) 1.179 (±0.049) 1.173 (±0.048) 1.134 (±0.030) 8.04 (±1.73) 7.72 (±1.52) 8.73 (±2.08) 9.65 (±2.17) 7.12 (±1.40) Table 21: Estimates of Kˆ well, Kˆ efu and `ˆ for the mean of all angular directions hh(r)i, single transects hh(r, j)) and five distinct combinations of heads with randomly chosen angular position hh˜ (r, Ck)ik=1,...,5. In brackets are given the 95% confidence intervals. The sampling size with regard to the radial distance remained constant at 1m intervals. hereafter referred to as transects hh(r, 1)i,..., hh(r, 4)) (Table 21). Then the drawdown curves were analyzed with a random combination of head measurements from any of the four axial directions hh˜ (r, C1)i,. .., hh˜ (r, C5)i, listed in Table 21. Secondly, we controlled the sample size regarding the radial distance to the pumping well, while using the average hydraulic head of all transects. The sample size was reduced with respect to equidistant sampling intervals for scenarios S2-S4 (Table 20). For scenarios S5 and S6 sampling intervals remained at 1m, but focused only on head measurements in the vicinity of the pumping well or the far field, respectively (notably, head measurement at the pumping well was excluded for S6). Likewise, head measurements focused on the vicinity of the pumping well or the far field for scenarios S7 and S8, respectively, but heads were sampled across the entire range of the aquifer (logarithmic scaled arrangement). Finally, we combined the two controlled experiments and randomly sampled head measurements from any of the four transects for all scenarios as defined in Table 20. 5.3.3. Parameter Estimation under Angular Limitation The drawdown curves in Figure 39 – angular mean, the four transects and one random angular combination of head measurements – were in good accordance in the far field and directly at the pumping well. Hence, estimates of Kˆ well and Kˆ efu were accurate and only little uncertainty was involved for both, the four transects hh(r, i)) and the random combinations hh(r, Ci)) (Table 21). The hydraulic conditions of both locations, pumping well and far field, were the same for all j and Ck, thereby resulting in similar estimates for Kˆ well and Kˆ efu. At a radial distance from the pumping well corresponding to one correlation length (i = 8m) the 5. Estimating Parameters of Aquifer Heterogeneity Using Pumping Tests Figure 39: Simulated drawdowns of a single realization (with statistics KG =10-4 m/s, 2 = 1, e =1, i =8m and h(R =32m)=0m). Black line denotes the angular mean hh(r)) of the four transects hh(r,1)i,...,hh(r,4)) (grey dashed lines). The black dots mark one random combination of head measurements hh(r,C1)) . drawdowns of the four transects hh(r,i)) (dashed lines in Figure 39) deviated from the angular mean hh(r)i, indicating that the estimated correlation length `ˆwas sensitive to the selected transect. Likewise, estimates of `ˆwere less accurate for selected transects, e.g., `ˆof hh(r,3)) = 12.52 in Table 21. A possible explanation for this is the occurrence of patches of much higher or lower than the average hydraulic conductivity in the environment of the pumping well, which critically influence hydraulic heads of the transition zone between Kwell-driven and Kefu-driven depression. For randomly combined head measurements from all four transects hh(r,Ci)) (black dots in Figure 39) the estimates of `ˆwere more accurate but less certain, as indicated by larger confidence intervals (Table 21). The random distribution of measurements around the pumping well buffers the effect of estimating local values as for single transects, but also results in a larger deviation from the mean head hh(r)) (Figure 39). Based on these findings, we conclude that for heterogeneous aquifers the angular arrangement of head measurements is critical to estimate sound and certain values of `ˆ. 5.3.4. Parameter Estimation under Radial Limitation For each scenario S1-S8 (Table 20) the estimates of Kˆ efu, Kˆ well and `ˆand the corresponding confidence intervals CI ˆ, CI ˆand CI`ˆ – all normalized by the initial parameters Kefu, Kwell KefuKwell and i – are depicted in Figure 40, where values close to one indicate high accuracy. Likewise, `˜ is plotted as alternative estimate to `ˆto depict estimates of the correlation length while Kˆ efu and Kˆ well being constant. For each scenario the approach of estimating `˜separately from Kˆ efu and Kˆ well was found to be very suitable to reduce the confidence interval CI`˜dramatically (Figure 40). 5.3. Pumping Tests under Limited Data Availability -a Paradigm for Field Applications 117 Figure 40: Normalized estimation results for eight scenarios of radially limited head data (Ta ble 20) with error bars marking the normalized 95% confidence intervals. Circles denote results for Kˆ efu/Kefu, squares for ˆ`/i and diamonds Kwell/hKwelli, stars for ˆ for ˜ `/`. Parameter estimates of Kˆ well and Kˆ efu were nearly insensitive to any sampling arrangement with regard to the radial distance except for S6 where the estimation of Kˆ well failed due to a lack of head measurement at the pumping well. As a consequence also estimates of correlation length were unreliable for S6. Otherwise, estimates of the correlation length were sensitive to increasing sampling intervals (S2-S4) and the distance of the head measurements from the pumping well (S5, S7 and S8). These findings emphasize that at least one (if measured at the pumping well) to three (if measured in the vicinity of the pumping well) head measurements are required to facilitate reliable estimates of Kˆ well. Likewise, four to five head measurements located in the far field of the depression cone enhance reliable estimates of Kˆ efu. Based on this findings, we conclude that limited sample size with regard to the radial distance from the pumping well affects the uncertainty of the estimates Kˆ efu and Kˆ well, represented by larger confidence intervals rather than their accuracy. This uncertainty also influences the uncertainty of the correlation length `ˆ, whereas the accuracy of `ˆ critically depends on the number of head measurements located at the transition zone (0 - 2`). This result is critical for field applications where the correlation length is the target rather than the input parameter. The combination of both approaches, sample size reduction with regard to angular and radial head limitations, even amplified the effects: the accuracy of parameter estimates decreased while confidence intervals increased, particularly for the correlation length. 5. Estimating Parameters of Aquifer Heterogeneity Using Pumping Tests Figure 41: Distribution of pumping wells (black dots) and observation wells (circles) in the area of investigation at the field site Horkheimer Insel (according to Schad [1997]). 5.4. Field Application Given our findings from the simulated pumping tests in section 5.3, we applied hefw(r) using empirical data from on-site pumping tests to estimate the parameters Kˆ G, ˆ2 and `ˆof a natural aquifer. 5.4.1. Field Site Horkheimer Insel We used empirical pumping test data from the Horkheimer Insel [Schad, 1997], located in the Neckar Valley in southern Germany. The aquifer consisted of unconsolidated fluvial sediments with poorly sorted sand and gravel and can be considered heterogeneous with a log-normal conductivity distribution [Schad and Teutsch, 1994]. On average the saturated thickness L was 3m, overlain by impermeable flood deposits and underlain by limestone formations. The dominating general hydraulic gradient was 0.001 towards the Neckar River. The infiltration rate from bedrock into the aquifer was negligible. According to the findings of Schad [1997] the average hydraulic conductivity derived from flowmeter measurements was KG =4.57e-3m/s. The average transmissivity was estimated as T =3.1e-2m2/s using Theis’ analytical solution, which is equivalent to Kefu = T/L =1.0e-2m/s. Detailed variogram analysis revealed that the vertical and horizontal correlation length ranged between `v =0.15 - 0.25m and f =6 - 10m, respectively, resulting in an anisotropy ratio of e = `v/f =0.015 - 0.04. The variance of hydraulic conductivity ranged from 2 =3.17 for permeameter measurements and 2 =2.32 for grain size analyzes to 2 =1.57 for flowmeter measurements. More details regarding the field site, wells, pumping tests, and parameter estimates of K(x) from laboratory trials are provided in Schad [1997]. For the purpose of this study we used data derived from small scale pumping tests (Figure 41). At each of the three pumping wells W40, W42, and W44, two pumping tests were performed. 5.4. Field Application For each of the pumping tests the drawdown was measured at four observation wells. The tests continued for two hours with pumping rates varying between 2l/s and 5.5l/s depending on the well yield. The pumping test duration of two hours is assumed to be ideal (neither too short nor too long) because (i) transient drawdown data at later times corresponds to a quasi steady state flow regime [Neuman et al., 2007], (ii) the influence of field-site dependent boundaries can be neglected until three hours of pumping [Schad and Teutsch, 1994], (iii) storativity, well bore storage and well loss have a significant influence to drawdown at early times which can be eliminated by using drawdown data at later times [Kruseman and de Ridder, 2000]; (iv) the influence of areas with low permeability do not affect effective conductivity until after two hours of pumping [Schad and Teutsch, 1994]; (v) after two hours the flow is essentially horizontal, whereas vertical flow becomes negligible and, hence, this phase is interpreted as influenced by inherent aquifer properties only [Schad and Teutsch, 1994]. 5.4.2. Analysis of Empirical Data The measured drawdown data hh(r)) of pumping tests PT40, PT42 and PT44 was compared to Thiem’s solution (Figure 42), such that it matched the drawdown near and at a large distance from the pumping well, eventually resulting in representative conductivities that can be considered first estimates of Kwell and Kefu, respectively. Then, hefw(r) was applied using nonlinear regression to estimate the parameters of K(x), as described in section 5.2.2. The boundary conditions at the pumping well were assumed to be BCH because pumping tests were performed at constant pumping rates and Kwell