Institute of Groundwater Studies, University of the Free State STOCHASTIC GROUNDWATER FLOW MODELS IN CONFINED AND LEAKY AQUIFERS Dissertation – MSc Degree Geohydrology Sarti Rautia Amakali 6/6/2019 Acknowledgement Foremost, I dedicate my utmost thankful to the Almighty God for grant me a chance again to work on this dissertation. Lord, you have been so much great to me all my life Iong, you blessed me with the right people, provided me with good knowledge and skills, and financial assistance, until the successful completion of this dissertation. My deepest appreciation goes to my supervisor, Professor Dr. AbdonAtangana, for his valuable guidance throughout the construction of this dissertation. Prof, you are so wonderful! I’m glad God grant me a chance to meet you and finally work with you on my dissertation. You pushed me to the right direction and I discovered a great confidence in myself, made me fearless of anything. Without your guidance and persistence Prof, this dissertation would not be completed by now. I would like to thank my family for their support. To my son, Junior, keep learning and focus boy, I need you to do more than what I delivered. You are always there for me. Not forgetting my friends, whom I regard as my biological family because they were so great of support in reflecting over my problems and findings, as well as providing happy distraction to clear my mind off my research, when I needed to rest and refresh my mind. I should thank the University of the Free State including the Post Graduate School for sponsoring me, allowing me to complete my studies here, regardless of the challenges I experienced which frustrated my efforts. 1 DECLARATION I hereby declare that the dissertation entitled STOCHASTIC GROUNDWATER FLOW MODELS IN CONFINED AND DETERMINISTIC MODELS IN CONFINED AND LEAKY AQUIFERS submitted by me in fulfillment of the Master‘s Science Degree majoring in Geohydrology from Institute of Groundwater Institute (IGS), is my original work and has not submitted ever before to Institute of Groundwater Institute or to any other institution for the fulfillment of any course of study. I also declare that no any chapter of this dissertation in whole or in part is incorporated from any earlier reports done by others or me. Name: Sarti Rautia Amakali Student no: 2015254739 Institute of Groundwater Studies University of the Free State 2 Table of Contents Acknowledgement .................................................................................................................................. 1 Abstract ................................................................................................................................................... 6 Keywords ................................................................................................................................................. 7 List of Symbols ........................................................................................................................................ 8 CHAPTER 1: INTRODUCTION ................................................................................................................... 9 1.1 Background of the study ......................................................................................................... 9 1.2 Problem statement ............................................................................................................... 11 1.3 Aims and objectives .............................................................................................................. 12 1.4 Research outline ................................................................................................................... 13 1.5 Research framework ............................................................................................................. 15 CHAPTER 2: CONFINED AND LEAKY AQUIFERS ..................................................................................... 16 2.1 Confined aquifer ......................................................................................................................... 17 2.1.1 Groundwater flow in confined aquifer ................................................................................ 20 2.2 Leaky aquifer ............................................................................................................................... 22 2.2.1 Groundwater flow equation in leaky aquifer ....................................................................... 27 CHAPTER 3: DETERMINISTIC AND STOCHASTIC MODELING ................................................................. 29 3.1 Deterministic models .................................................................................................................. 31 3.2 Stochastic models ....................................................................................................................... 32 3.3 Literature review on applications of stochastic and deterministic models ................................ 33 CHAPTER 4: DISTRIBUTION .................................................................................................................. 39 4.1 Statistical data analysis ............................................................................................................... 40 4.1.1 Mean .................................................................................................................................... 40 4.1.2 Variance ............................................................................................................................... 41 4.1.3 Standard deviation ............................................................................................................... 41 4.1.4 Skewness .............................................................................................................................. 42 4.1.5 Kurtosis ................................................................................................................................ 42 4.2 Density distributions ................................................................................................................... 42 4.2.1 Normal /Gaussian distribution ............................................................................................. 43 4.2.2 Binomial distribution ............................................................................................................ 44 4.2.3 Lognormal distribution......................................................................................................... 44 4.2.4 Binomial–Exponential 2 distribution .................................................................................... 45 4.2.5 Poisson distribution ............................................................................................................. 45 4.2.6 Geometric distribution ......................................................................................................... 45 4.2.7 Exponential distribution ....................................................................................................... 46 4.2.8 Generalized exponential distribution .................................................................................. 46 3 4.2.9 Bernouili distribution ........................................................................................................... 46 4.2.10 Weibull distribution ........................................................................................................... 47 4.2.11 Beta distribution ................................................................................................................ 47 4.2.12 Gamma distribution ........................................................................................................... 47 4.2.13 Rayleigh distribution .......................................................................................................... 48 4.1.14 Cauchy distribution ............................................................................................................ 48 4.1.15 Chi‐squared distribution .................................................................................................... 48 4.1.16 Empirical distribution ......................................................................................................... 48 4.1.18 Pareto distribution ............................................................................................................. 48 4.1.19 Mittag‐leffler distribution .................................................................................................. 48 4.3 Estimation methods .................................................................................................................... 49 4.4 Hydraulic parameters .................................................................................................................. 50 4.4.1 Transmissivity ....................................................................................................................... 51 4.4.2 Storativity ............................................................................................................................. 53 4.4.3 Leakage factor ...................................................................................................................... 54 CHAPTER 5: ANALYSIS OF STOCHASTIC MODELS OF GROUNDWATER FLOW: CONFINED AND LEAKY AQUIFERS .............................................................................................................................................. 56 5.1 Analysis of stochastic model of groundwater flow: Confined aquifers ...................................... 57 5.2 Analysis of stochastic model of groundwater flow: Leaky aquifers............................................ 59 CHAPTER 6: APPLICATION OF NEWTON METHOD ON STOCHASTIC GROUNDWATER FLOW MODELS FOR CONFINED AND LEAKY AQUIFERS ................................................................................................. 61 6.1 Application of Newton method on stochastic Theis’s confined aquifer ..................................... 66 6.2: Application of Newton method on stochastic Hantush’s leaky aquifer .................................... 68 6.3 Stability for the stochastic confined aquifer equation ............................................................... 71 6.4 Stability for the stochastic leaky aquifer equation ..................................................................... 79 6.5 Simulation ................................................................................................................................... 86 CONCLUSION ......................................................................................................................................... 98 REFERENCES ........................................................................................................................................ 100 4 List of figures Figure 1: Illustrate of research framework covered in this dissertation……..…………………………15 Figure 6.1: Numerical solution with harmonic mean for transmissivity and storativity ...................... 87 Figure 6.2: Contour plot for harmonic mean ........................................................................................ 88 Figure 6.3: Numerical simulation with arithmetic mean ...................................................................... 89 Figure 6.4: Contour plot with arithmetic mean .................................................................................... 89 Figure 6.5: Numerical solution with geometric .................................................................................... 90 Figure 6.6: Contour plot of geometric mean ........................................................................................ 90 Figure 6.7: Numerical solution with geometric mean .......................................................................... 91 Figure 6.8: Contour plot with geometric mean .................................................................................... 91 Figure 6.9: Numerical simulation with arithmetic mean ...................................................................... 92 Figure 6.10: Contour plot with arithmetic mean .................................................................................. 92 Figure 6.11: Numerical simulation for harmonic mean ........................................................................ 93 Figure 6.12: Contour plot with harmonic mean ................................................................................... 93 Figure 6.13: Numerical simulation with harmonic mean ..................................................................... 94 Figure 6.14: Contour plot with harmonic mean ................................................................................... 94 Figure 6.15: Numerical simulation for arithmetic mean....................................................................... 95 Figure 6.16: Contour plot for arithmetic mean .................................................................................... 95 Figure 6.17: Numerical simulation for geometric mean ....................................................................... 96 Figure 6.18: Contour plot with geometric mean .................................................................................. 96 5 Abstract This dissertation proposes an application of stochastic modeling of groundwater flow in confined and leaky aquifers. We are estimating that aquifer parameters such as transmissivity, storativity and leakage factor vary, not constant, in space at different period especially in heterogeneous environment. Heterogeneous environment are known to be complex because of their uncertainty. Uncertainty referred in modeling includes errors in dataset, which might be bias or variance (under fitting/over fitting), low or not enough data, or unbalanced data, which all affect the model produced if not captured with appropriate model technique. The groundwater flow equation for confined and leaky aquifers derived by the latest version Atangana and Ramotsho, as well asAtangana and Mathobo, which all include scaling matrix of the soil, are considered and further modified to a new scheme of stochastic models for confined and leaky aquifer. We tried to achieve the capture statistical setting of aquifer parameters using the concept of stochastic modeling technique. The aquifer parameters are replaced by distribution for instance, Gaussian or normal distribution.Due to the complexity of the modified models, it is almost impossible to obtain the exact solution by using analytical solution, thus we opt to numerical analysis, in particular the Newton method used to derive the numerical solutions of the modified models. Detailed analysis of stability and convergences, we used method are presented for both models. Numerical simulations are depicted for different distributions. 6 Keywords Arithmetic mean Complex environment Confined Aquifer Crank Nicholson Deterministic Discretization Distribution First Order Derivative Geometric mean Groundwater flow model Harmonic mean Heterogeneity Leakage Leaky aquifer Newton method Numerical simulation Second Order Derivative Stability Statistic Steady state flow Stochastic model Storativity Transmissivity Unsteady state flow 7 List of Symbols BC: Boundary Conditions CDF: Cumulative Density Function DK: Delaunay Kriging DSHW: Double Seasonal Holt Winters EKF: Extended Kalman Filter FD: Finite Difference FE: Finite Element GP: Gaussian Processes GLUE: Generalized Likelihood Uncertainty Estimation IDW: Inverse Distance Weighting LU: Lower Upper ODE: Ordinary Differential Equation OK: Ordinary Kriging MC: Minimum Curvature PDE: Partial Differential Equation PDF: Probability Density Function REV: Representative Elemental Volume SSLE: Sequential Successive Linear Estimator SG: Sequential Gaussian UK: Universal Kriging 8 CHAPTER 1: INTRODUCTION 1.1 Background of the study Effective groundwater management is becoming a compulsory tool in almost every industry due to accelerating demand of groundwater resources resulting in its depletion and limitations. Groundwater has become extremely important, pushing governmental institutions and non‐governmental institutions to work together to set up legislations to manage groundwater resources sustainably. In the past, most industries relied on surface water due to its availability that time, but nowadays most surface water sources such as rivers, streams, canals, and springs are heavily polluted and depleted. This has pushed much demand on groundwater use. Driest countries like Namibia receives little rainfall annually (±400mm) to feed surface water basins such as rivers, streams and pans (Kumba, 2003), thus relies much on groundwater use, which needs to be well managed. The strategies used by water service providers in allocating groundwater quotas to the farmers and other bulk water users should consider the long‐term aquifer sustainability and potential effects of abstraction on the reliability of water supply. Groundwater occurs in the subsurface, and we cannot easily estimate its quantity and quality because we unable to see it physically. Scientists have come up with different modeling techniques to estimate hydraulic parameters in an aquifer, which helps to determine the occurrence of groundwater in a system. The parameters are usually obtained from aquifer pump tests data, which yield information about the specific aquifer’s transmissivity, storativity, leakage or hydraulic conductivity. The geological properties of the site also need to be identified as they can help the modeler to decide which applicable modeling technique will be suitable for such system. Some aquifers have 9 uniform geologic stratification which needs to be treated as homogenous, while some have complex geologic structures need to be treated as heterogeneous in groundwater modeling. But the challenge is often observed indifferent sedimentary basins with their geologic formations vary in space and time, which need to be treated as heterogeneity in modeling. Analyzing a heterogeneity environment is quite a challenge in terms of groundwater flow because heterogeneity occurs at all scales of the formation, and needs to be investigated with appropriate modeling technique. There are different approaches of modeling widely used by scientists to address problems to fit their needs. Empirical black box, probabilistic, deterministic and stochastic models are common model types used nowadays in groundwater flows to describe the flow behavior of groundwater systems (Adem and Batelaan, 2006). Empirical black box models are obtained from experimental data supported with mathematical functions like law of Darcy. They are known to be good in simulating spatial variability of the subsurface. Black box models deliver a mathematical mapping between historical inputs and outputs without requiring physical information on the investigated system. In water sector, black box models are suitable for nonlinear, non‐Gaussian and non‐stationary processes relate to the research. Probabilistic models use statistics and elements which also accepted by stochastic models. One can start with a simple probabilistic distribution model type, but may end up with a complicated stochastic time dependent model. Deterministic models are based on fixed mathematical descriptions of natural processes. They are commonly used in groundwater flow investigations as they are easy to use and can predict a single outcome from a given set of 10 circumstances. In short, they represent the physical processes we observe in the real world resembling the subsurface flow. For deterministic modeling, model parameters should be at least known and uniquely determined in the model by its variables. Stochastic models are based on mathematical description of natural processes, but include statistical analysis (probabilistic). Generally, stochastic is originated from a Greek term σᴛόxo (stόkhos) which means guess (Renard, 2007). Stochastic models are described by random variables, which are usually guessed. In this dissertation, we use stochastic modeling to estimate transmissivity, storativity and leakage, which vary in confined and leaky aquifers. The techniques used to solve governing equations in modeling include Finite Difference (FD) and Finite Element (FE). We are applying Finite Difference in our work to solve numerically Partial Difference Equation (PDE) to obtain an exact solution, which will be used in our stochastic models for confined and leaky aquifers. The numerical model will provide approximate solutions to the governing equation via discretization of spatial (space) and temporal (time). We noticed that hydraulic parameters are not constant in space at a time; they keep changing, so we are treating them as distribution probabilities in a form of normal, log normal distribution etc. to find the exact form of distribution fitting each parameter. 1.2 Problem statement On this dissertation, we used stochastic modeling to investigate hydraulic parameters such as transmissivity, storativity and leakage factor in confined and leaky aquifers. The strength of stochastic models that caught our interest is that we are able to address uncertainty in heterogeneous environment. 11 Uncertainties may arise from a lack of knowledge of the actual distribution of subsurface system, under or overestimation of discharge rate, recharge rates or hydraulic heads assigned during the calibration stage. It depends on the choice of the modeler to select the type or kind of model to use, but s/he must ensure that the model preferred should be well designed and adequately resemble the natural system it represents (Fetter, 2001). Groundwater flow models use elements described by equations representing the hydraulic parameters including the physics, properties, stresses and geometry of the aquifer system. The mathematical equations solved will provide the hydraulic parameters corresponding to the conceptual model used to simulate the aquifer. In this dissertation, a numeric type of stochastic approach with the use of partial differential equations will be used to capture the movement of subsurface water within confined and leaky aquifers. Most of the past groundwater modeling studies focused on spatial distribution of hydraulic conductivity, had used deterministic models which can be solved analytically or numerically, by using partial difference equations with finite differences methods (Dagan, 1989; Gelhar, 1993; Hsu, Zhang and Neuman, 1996). 1.3 Aims and objectives The aim of this research is to develop a stochastic model of groundwater flow in confined and leaky aquifers. The objectives: a) To capture statistical setting of aquifer parameters using the concept of stochastic modeling technique. 12 b) To modify the mathematical equation describing the flow of subsurface water within a confined and leaky aquifers by converting the constant aquifer parameters to probability distributions. c) To analyze numerically the modified models using a new scheme (Newton method). d) Use of stochastic models compared to other simulation methods such as deterministic models in groundwater flow systems. 1.4 Research outline The dissertation has six chapters. Chapter 1 outline the introduction of the dissertation including the background of the proposed study, problem statement explaining the importance of such research needed to be carried out. The aims and objectives of the dissertation are also included in chapter 1. The following chapter 2 focused on confined and leaky aquifer in general. The Theis (1935) equation for unsteady state conditions and Thiem (1906) equation for steady state conditions in confined aquifers is included in this chapter. The mathematical equation for leaky aquifer developed by Hantush & Jacob (1955) is also included in this chapter. The existing model for groundwater flow in confined aquifer as derived by Theis (1935) is discussed and modified to a new equation we developed which included the soil matrix which was neglected in the previous model. The mathematical model for leaky aquifers is also described in this chapter. In Chapter 3, we discussed stochastic and deterministic models in general including literature review including their applications on groundwater flow within subsurface. Chapter 4 covers distribution probabilistic, explaining how hydraulic parameters such as transmissivity, storativity and leakage can be obtained through means of statistics and distribution methods. Chapter 5 entails the analysis of a new 13 stochastic model in confined and leaky aquifer. In Chapter 6, we introduce the new scheme of Newton method being applied to stochastic groundwater flow equation of confined and leaky aquifers obtained in the previous chapter 5. It also entails the conversion of 1st order derivative and 2nd order derivative and the concept of discretization for Partial Difference Equations into the newly formed equation. Stability for all numerical equations of confined and leaky aquifers is also discussed in this chapter. Thereafter, simulation for the models concluded this chapter. 14 1.5 Research framework The research framework used on this dissertation to achieve the aims and objectives is described in figure 1 below: Reviewing literature on deterministic and stochastic models of groundwater flow in confined and leaky aquifers Use of distribution probabilistic and statistical analyses for hydraulic parameters, which varies in space at a time De rivation and analysis of a new Derivation and analysis of a new stochastic model for confined aquifers stochastic model for leaky aquifers Derivation of Stochastic model in Derivation of Stochastic model in leaky conf ined aquifers aquifers Application of Newton method on Application of Newton method on stochastic Confined model stochastic leaky model Stability for the confined model Stability for the leaky model Simulation for the confined model Simulation for the leaky model Figure 1: Illustrate of research framework covered in this dissertation 15 CHAPTER 2: CONFINED AND LEAKY AQUIFERS Groundwater is found beneath the earth surface, stored within the rock pores and fractures. When it rains, some part of the surface water infiltrate through soil and percolate into underground, becoming groundwater resource. Streams, springs and rivers can lose water (effluent) to the underground water to recharge the aquifer. Groundwater can stay safely beneath earth surface in units called aquifers for years and decades. Aquifers are grouped as confined, unconfined and leaky aquifers. In this dissertation, we only focus on confined and leaky aquifers. Generally, the aquifer can be characterized by the nature and distribution of its lithology, flow behavior and geological structures and their formations (Payne and Woessner, 2010). Lithology, which obtained through drilling of wells, describes the aquifer’s physical logs that consist of mineral composition, grain size, sorting and packing of geological features describing an aquifer system. Stratigraphy characterizes the geometry and age of formations of geological system. Geological structures refer to secondary features such as fractures, faults, folds that may occur after deposition. The response of aquifer is controlled by the hydraulic parameters namely permeability, hydraulic conductivity, transmissivity, storativity and leakage factor. The fluid materials (water) determine permeability, its ease for the water to move through the aquifer materials. Hydraulic conductivity is the ability of the aquifer materials to let water to pass through it. Transmissivity, storativity and leakage factor are the parameters of interest in this dissertation, and will be discussed in the following sections. 16 2.1 Confined aquifer In Fetter (2001), a confined aquifer is described as an aquifer that is layered below and above by impermeable rocks (aquitard) which limits the groundwater movement into or out of the aquifer. It means water found on this aquifer types is bounded within the aquifer, only leaves the aquifer by means of discharge through drilled wells. Confined aquifers are usually found deep underground. Groundwater in confined aquifers is often found to be under pressure, rising to the potentiometric level because of a layered aquitard. The aquitard is beneficial as it protects the aquifer from contamination because polluted water cannot easily enter the confined aquifer. The aquifer gets recharge at points where it is exposed to the surface by means of rain or loose streams, rivers and springs that cross the outcrops and deep underground tributaries. Water moves slowly along the dip of the strata below the overlying impermeable layer to the aquifer. The cone of depression caused by pumping in confined aquifers can be described as shallow but widely, extending at a fast speed and fill up at a slow rate. Thiem’s method (1906) (steady‐state flow) and Theis’s method (1935) (unsteady) are the solutions used to determine aquifer parameters in confined aquifers. According to Kruseman, De Ridder, and Verweij (1970), the assumptions of confined aquifers for steady state flow (Thiem) conditions: 1. At the beginning of pumping, the potentiometric surface is constant over the area affected by pumping; 2. The aquifer is confined at the top and bottom with apparent infinite extent; 17 3. The aquifer is homogenous, isotropic, and of a uniform thickness over the area influenced by pumping; 4. The well is pumped at a constant discharge rate; 5. The well fully penetrates the aquifer’s thickness and water flows horizontally; 6. The aquifer recharge through areas where it’s exposed to the surface and via leakage from other layers. Unsteady-state flow conditions (Theis, 1935): 7. The aquifer can compress, and water is discharged immediately from the storage due to a drop in head; 8. The well diameter is small so that well storage is negligible. Pumping in confined aquifers affects a large area (distance) of the aquifer, thus its cone of depression always extend far further away. Therefore, it is always advisable to place observation wells or piezometric tubes at reasonable distances from the pumped well to measure drawdown. In case of absence of observation well, for analytical analysis one has to accommodate an imaginary well together with a real well to get the hydraulic head correctly. The observed head will help to calculate aquifer parameters with the application of analytical solution or numerical solution which can be Theis (1935) or Thiem (1906) (Hu, Chen and Chen, 2011). The mathematical equation (1) commonly used for confined aquifer as derived by Theis (1935) for calculating hydraulic head in unsteady state flow conditions is as follows: 18 Q − (1) ℎ − ℎ = du 4πT u Integral can be replaced in this way as indicated in equation (2) below: (2) (,)= () 4 where; h0: initial head h: head s(r,t): drawdown at distance (r) at time (t) after the start of pumping (L) Q: discharge rate (L/s), it shall be converted to m3 per day even the well is pumped for less than 24 hours T: transmissivity W(u): well function uobtained from equation (2) above can be solved with this formula in equation (3) below: ² (3) = 4( − ) where; r: radial distance from pumping well (m) S: aquifer storativity (dimensionless) t: elapsed time since pumping began 19 t0: initial time pumping T: transmissivity Thiem (1906) applied equation (4) for steady state flow conditions for calculating hydraulic head: 2 ( − ) (4) = where; Q: well discharge KD: transmissivity of aquifer r1, r2: respective distances of the piezometers from pumping well (m) s1, s2: respective steady‐state drawdown in the piezometers (m) 2.1.1 Groundwater flow in confined aquifer Groundwater flow equations are used to represent the properties of a system as they are assumed always constant. This is what is called Representative Elemental Volume (REV). In this dissertation, we are challenging this concept as the aquifer properties may not always constant in space at a time due to variability in spatial distribution which can relate to the inconsistent of geometry or pore sizes of aquifer system. There can be changes in aquifer’s geological structures, hydrologic cycles including precipitation to recharge, stream flow, and human induced activities such as over abstraction at a particular point of an aquifer. Nevertheless, most authors/modelers ignored this, so we are arguing that the aquifer parameters cannot be modeled as constant effectively. Also the drawdown used may only affect the pumped well and observation wells considered which only represent several points (part) of 20 the aquifer, not responding to the whole aquifer system, unless it’s a homogenous environment as most hydrologists and modelers concluded in their findings. As derived from Theis (1935), the groundwater flow equation (5) in confined aquifers described below determines the movement of groundwater in a radial distance of a system (space) at a particular time: (5) ℎ(,) ℎ(,) 1ℎ(,) = + where, S represents aquifer storativity (dimensionless), T for transmissivity (m²/d) in a radial distance (r), and (t) is the pumping time. The above equation does not take into account the structure of the geological formations; it considers only the transmissivity and the storativity of the aquifer, which in the normal situation cannot be representative as the water move slowly within the aquifer passing through each portion of the matrix soil. Thus neglecting the scaling effect will probably lead to misleading results. To solve this problem, Mathobo and Atangana (2018) suggested a revision of the Theis model starting from the Darcy´s law and using the mass balance equation, without neglecting any terms as was done by Theis. They suggested a modified equation (6), which is apparently more complex and more informative than the existing model. Their model includes the scaling effect of the matrix soil and is given as: ℎ(,) ℎ(,) ∆ ℎ(,) (6) = 1 + + 21 where, S represents aquifer storativity (dimensionless), T for transmissivity (m²/d) in a radial distance (r) at a time (t). In our work, we will consider the above equation and transform it to a stochastic model for further investigation. 2.2 Leaky aquifer Leaky aquifers, known as, semi confined aquifers are those aquifers types partly bounded by aquitard (semi permeable) on its upper and lower boundaries, or in some cases the upper boundary can be aquitard while the lower boundary can be aquiclude. ŞEN (2000) also agreed that leaky aquifers are much complicated to investigate because additional water from upper aquitard may enter the aquifer any time nonstop. So modeling this aquifer type, one has to be vigilant with the effects of leakage factor. The geological materials that made up the overlying or underlying aquitard (confining layer) can be impervious but can compress/bend leading to a significant leakage through the layers, which in some cases ignored. As derived from Kruseman, de Rider and Verweij (1991), the following are assumptions of leaky aquifers in steady state flow conditions: 1. The aquifer is leaky and has an seemingly infinite areal extent; 2. The aquifer is bounded on top by aquitard; 3. The aquitard is overlain by an unconfined aquifer called source bed; 4. The water table in source bed is initially horizontal; 5. The water table in the upper aquifer (source bed) does not drop during pumping of the aquifer as it receives water from outside source; 22 6. As pumping continues, water table in the upper aquifer will fall, as more of its water will be leaking through the aquitard into the pumped aquifer. 7. Groundwater flow (leakage) in the aquitard is vertical; 8. The aquitard is incompressible, so that no water is released from storage in the aquitard when the aquifer is pumped; 9. The aquifer and the confining layer (aquitard) are homogenous, isotropic, and of uniform thickness over the area influenced by pumping; 10. Prior to pumping, piezometric surface and water table remains horizontal over the area influenced by pumping; 11. The drawdown in the unpumped aquifer or in the aquitard is negligible. 12. The well is pumped at a constant rate; 13. The well fully penetrate the aquifer thickness and receive water by horizontal flow; Unsteady-state flow conditions: 14. The aquifer is compressible, and water released from storage in the aquifer and the water supplied by leakage from aquitard is discharged instantaneously with a decline in head; 15. The well diameter is small so that well storage is negligible. The Hantush and Jacob (1955)developed mathematical analysis for leaky aquifer used to solve problems for fully penetrating well in leaky confined aquifer of infinite extent using the following equation (7): 23 / (7) = 4 Integral can be replaced and written as in equation (8) below: (8) = where, b’: aquitard thickness (L) K’: vertical hydraulic conductivity of the aquitard (L/T) Q: pumping rate (L3/T) r: radical distance from pumping well to observation well (L) s: drawdown t: elapsed time since start of pumping (T) T: transmissivity (L2/T) S: storativity When pumping from a leaky aquifer, at the beginning the groundwater flow is assumed horizontal but as pumping goes on, the flow becomes vertical in underlying and overlying aquitard. When the flow becomes vertical, the well is now supplied by storage from aquitard. The drawdown considered here is a product of pumping rate, transmissivity and well function. The drawdown in a leaky aquifer is assumed proportional to the leakage rate from aquitard over the entire period of pumping. If there is sufficient leakage resulting from 24 overlying or underlying aquitard of a pumped aquifer, the drawdown will be minimal. The size of the aquitard whether it is thin or thick does matter as long it affect the aquifer’s compressibility. Thin aquitard accompanied with long pumping tests may result in low storativity when calculated with Hantush Jacob formula. Thick aquitard and short pumping durations may result in high storativity when Hantush formula is being used. Hantush and Jacob (1955) developed a solution to analyze the condition of thick aquitard using the equation (9) below: (9) But this solution (9) did not take into account the storativity effects, and then later it was modified by Hantush (1960) to include the storativity ending up the solution gives zero as shown in equation (10) below. The solution has been criticized so far in leaky system as it is not always clear whether leakage comes from above the aquitard of the pumped aquifer or below aquitard. It takes to consideration that if there is only significant leakage or pump test did not last for long, the results will conclude: = 0 (10) It means that there is no leakage observed, but in reality there is significant leakage and it is usually ignored. For a better analysis of leakage factor, Neuman and Witherspoon (1969) developed a solution method, which needs an observation wells located in the aquitard below and above, and in the aquifer being pumped. The observed drawdown in the aquitard will be used in a ratio to that drawdown measured 25 in the pumped aquifer at the same time and same radial distance from pumping well to estimate hydraulic properties of an aquitard. For partially penetrating observation well in a leaky confined aquifer, the following equation (11) applies, as derived from Hantush (1960): (11) = (,/) 4 2 + (− ) − . . (,(/)²+ / ² where; b: aquifer thickness (L) d: depth to the top of pumping well screen (L) d’: depth to the top of observation well screen (L) Kr: radical (horizontal) hydraulic conductivity (L/T) Kz: vertical hydraulic conductivity (L/T) l: depth to the bottom of pumping well screen (L) l’: depth to the bottom of observation well screen (L) (u,B): the Hantush and Jacob well function for leaky confined aquifers (dimensionless) z: piezometer depth (L) 26 2.2.1 Groundwater flow equation in leaky aquifer In the last decade the groundwater flow model within a leaky aquifer was suggested by Hantush and Jacob (1955) used intensively to model such physical problem. The equation is derived based on the Theis model, which is also borrowed from the heat equation. The groundwater flow equation (12) in leaky aquifer derived from Hantush and Jacob (1955) model, when specific storage of confining unit is assumed negligible and the head in overlying aquifer is assumed unaffected: ℎ 1 ℎ ℎ ℎ (12) + − = Hantush (1960) modified the equation (12) to include the specific storage of confining unit and derive the following equation (13): ℎ 1 ℎ 1 ℎ(,) ℎ + − = (13) where, B: dimensionless leakage parameter h: drawdown r: radial distance t: time z: vertical distance in confining layer Although the model has been used with great success in the last past years, one need to note that, some assumptions or modification and simplifications were done to obtain a simple model. Amanda and Atangana (2018) undertook an investigation of the groundwater flowing within a leaky aquifer; they reused the Darcy law and the mass balance concept to revisit the existing model. 27 Without neglecting any parameters and simplification of any terms, they produced a new mathematical model (14) if groundwater flowing within a leaky aquifer. The obtained mathematical model is given below as: ℎ ℎ ∆ 1ℎ ℎ(,) (14) = + + + ƛ² where, S represent aquifer storativity (dimensionless), T for transmissivity (m²/d) in a radial distance at a time. The above equation has a new parameter, which takes into account the scaling of the geological formation, the model does not only depend on the transmissivity, storativity, leakage factor but it depends now on the scale of the geological formation. In this section we use the model suggested by Amanda and Atangana (2018), and further transform it to stochastic equation to be able to capture above the leaky and the scale factor, the heterogeneity associate to statistical setting of the aquifer parameters. 28 CHAPTER 3: DETERMINISTIC AND STOCHASTIC MODELING Stochastic and deterministic approaches are modeling techniques that can be used to simplify physical structures of a real system about the quantity and quality of groundwater. These types of models help to understand the behavior of groundwater flow within the subsurface. The model requires data serving as input functions, which can be obtained from field pumping test data with the aim of getting hydraulic head to estimate transmissivity, hydraulic conductivity, storativity, leakage factor etc. The main difference between these approaches is the ability to handle errors in a dataset. The errors can be introduced in a dataset by means under fitting or over fitting, not enough or low quality data, bias and variance, and unbalanced data. All these error type will lead to poor quality model produced with not enough features to represent the actual system. Most of these errors arise from heterogeneous environment. Deterministic model cannot handle most of these errors, thus they are always applied to homogenous environment. Stochastic has the ability to handle most of the errors even in complex heterogeneous environment. Deterministic is a commonly used approach in most homogenous environments as it is simple to use and can predict a single outcome from a set of data, because it uses known variables. Its weakness lies on its limitations, as it is unable to address uncertainty in heterogeneous environment. Deterministic approaches have been used to model different scenarios because they are regarded as simple to use and not so complex when compared to stochastic. Groundwater models obtained in deterministic approach do not provide unique solutions as their combinations of variables may produce similar results (Renard, 2007).Therefore, stochastic models are more recommended in 29 quantitative research to run groundwater flow models because they are unique and their strength to take random values to solve errors in heterogeneous environments. Stochastic predicts multiple set of possible outcomes weighted by its likelihood or probabilities. In addition, the algorithms produced by stochastic models are independent representing random functions (van Leeuwen et al., 1998). The problem with stochastic modeling is that they are more complicated to use as they require a random number which needs to be obtained via simulation methods such as Monte Carlo, variance, kriging, and quasi‐Newton. According to modelers, they found out that it is difficult to use and time consuming. In terms of qualitative, studies reviewed by Wang, Ocampo‐Martinez and Puig (2016) shows better quality models obtained through stochastic modeling than those obtained in deterministic. The strength of stochastic groundwater flow models is that they are more realistic and non‐uniform. Deterministic and stochastic approaches use different interpolation methods depending on the project type and its needs. For deterministic, Inverse Distance Weighting (IDW), Minimum Curvature (MC) and Radial Basis Function (RBF) methods are used for interpolation. Interpolation methods used in stochastic includes Ordinary Kriging (OK), Universal Kriging (UK), and Delaunay triangulation kriging (DK) (Adhikary and Dash, 2017).Ordinary kriging can be used in fixed data sets to estimate the values of points. OK assumes a random function with a constant unknown mean value, and their estimate depends on semivariogram. Universal kriging is used for irregular data taken at various points in space resulting in trends. UK uses algorithm estimation results, and requires the drift function and semivariogram of the residuals. DK uses Delaunay triangles to determine search neighborhood around the estimation point. Semivariogram works well with isotropical model like exponential, 30 gaussian, spherical, power law, linear and Spartan semivariogram models (Varouchakis and Hristopulos, 2013). When estimating values, for example transmissivity, at unsampled locations, one can also use both measured transmissivity and hydraulic head values or specific data capacity instead of only using transmissivity. Specific capacity and hydraulic head are always correlated, thus can be used together. 3.1 Deterministic models Deterministic models are defined by differential equations with known number resulting in known output for well‐defined linear models with multiple outputs possible for non‐linear outputs as well. Equations are solved by numerical methods such as Finite Difference. Constructing a model through a deterministic technique proves to predict about what to expect in future and the influence of certain parameters in an aquifer system (Renard, 2007). Certain parameters such as transmissivity or storativity can vary at a point due to the geometry or geological structure of an aquifer as discussed on the previous chapter. The prediction about aquifer parameters from a (numerical model) is very necessary for hydrologists to understand groundwater flow within the aquifer. Koch and Arlai (2007) also agree that in deterministic models, the parameters used in the equation to represent the system are constant and even measured at several locations in aquifer system, resulting in errors which is usually goes unnoticed, and it is a challenge in heterogeneous environment. For instance, concern raised in areas where data measurements were not completed as hydrologists treated the system as continuum and conclude that it is a homogenous environment, which is not always true, might be heterogeneous and they modeled it as deterministic. We should always 31 remember that a model found to be inaccurate will mislead the audience and decision makers and affect groundwater resource. 3.2 Stochastic models Errors in a dataset results from poor data collection methods, unknown boundary barriers or geometry which all serve as input for the model and influence the quality of the model produced, if not well represented. This is exactly what affects deterministic models due to its inability to cover these types of errors. Wang, Ocampo‐Martinez and Puig (2016) stated that stochastic models are capable of combining all forms of uncertainty influencing the data by replacing them with random values representing the hydraulic parameters. Stochastic approach is a solution to model complicated environment, whereby the available data has to be averaged by means of statistical analysis and distribution probabilities, like what we are doing on this dissertation, to estimate effective aquifer properties that varies in an aquifer system. Heterogeneity environments are well represented by stochastic models because any uncertainty arising in the aquifer system can be captured. As like for deterministic models, stochastic processes can also assume the variation of a system in future, but most importantly, it allows the modeler to test the validity of the assumptions through statistical analysis and estimate the expected value of the future as well as the variation of that expected value (Sinha, and Prasad, 1979). It also allows one to monitor the hydraulic parameters predicted in a model within the range one would expect. 32 3.3 Literature review on applications of stochastic and deterministic models Different scenarios relate to the applications of stochastic models compared to deterministic are discussed below. Use of stochastic technique with Bayesian Inference Stochastic models use one or more elements, which can be solved analytically or numerical. Stochastic model based on Gaussian Processes (GP) addressed a fault tolerant problem using a Bayesian Inference theory when it was applied to the Barcelona drinking water network (Wang, Ocampo‐Martinez and Puig 2016). Bayesian Inference is particularly useful in stochastic modeling technique especially when not enough data available and required to deliver accurate results in every area of predictive system. GP regression algorithms were combined with Double Seasonal Holt Winters (DSHW) to generate disturbance forecast uncertainty into the system in space and it yielded good results. The simulation for the model worked well, managed to handle infeasibility occurred in the system. Deterministic and Stochastic models on groundwater saline contamination In Bangkok, Koch and Arlai (2007) investigated a heavily stressed saline contamination of a multi aquifer system. They used both deterministic and stochastic techniques to find out which model provide better accurate results representing the real system. The multi aquifer consists of a topmost soft clay layer and eight water bearing layers, and it recharges at the basin flanks where 2nd and 4th aquifer outcrops. Due to groundwater over abstraction on the 2nd, 3rd and 4th aquifers, the water levels declined so badly in some areas creating a cone of depression, leading to land subsidence and salt water encroachment in 33 some places due to sea water intrusion or vertical seepage of saline connate from marine clay layer. They used pumping and piezometric head data to construct the model. When they used deterministic approach, piezometric head data could not be represented perfectly, leaving a non‐zero residual. This error occurred because piezometric and pumping data are not exactly measured and deterministic calibration parameters obtained only represent local instead of a global minimum piezometric surface response. This type of aquifer can be assumed to be heterogeneous, and it can only represented well by stochastic approach. When stochastic approach was performed by means of assigning random field data using logarithmic transmissivity Y=lnT field with various sets of variances σγ2 and correlation lengths (,) for each layer which defines stochastic range in a multi‐layer aquifer system simulated. Stochastic modeling used Monte Carlo Modflow simulation to investigate how variance σ 2Y contributed to σ 2 H of observed head and residuals. The investigation results indicates that σ 2H and σ 2 Y are related to each other asσ 2‐ H σ 2 2 Y * λ , just as stochastic theory predicted. So in this scenario, stochastic model yielded accurate results than the deterministic approach does. Stochastic model of slip spatial amplitude of the fault Another scenario of stochastic modelling, done by Lavalle´e and Archuleta (2003) investigated a heterogeneous spatial distribution of the slip or pre ‐ stress along the fault zone of the 1979 Imperial Valley earthquake in California, USA. They applied stochastic modeling based on Gaussian distributions and non‐Gaussian distributions. They estimated power spectrum (as random variables) and interpolate the slip distribution. Probability Density Function (PDF) associated with power spectrum was also estimated. Fourier transform 34 and Inverse were used to represent the slip amplitude and position of the strike. It resulted that stochastic model based on non‐Gaussian distributions described well the spatial variability of the slip amplitude over the fault. The stochastic model based on Gaussian distribution failed to reproduce the spatial variability observed in the original slip distribution. The results proved that non Gaussian distribution works better with stochastic models. Groundwater level in a sparsely monitored basin The study conducted by Varouchakis and Hristopulos (2013) compares the performance of stochastic and deterministic methods for mapping groundwater levels in a sparsely monitored basin. For deterministic, they used Inverse Distance Weighting (IDW) and Minimum Curvature (MC) as interpolation methods for groundwater level mapping. IDW is an exact, fast, straight forward and computationally non‐intensive. Its limitation is that it cannot handle uncertainty measure, and arbitrary choice of the weighting function. MC is based on reduction of the total square curvature of the surface subject to the data constraints. Its weakness is that it suffers from oscillations caused by outliers in the data or very large gradients, especially in cases of small data. Deterministic methods use closed form mathematical formulas for IDW or solution of linear system of equations to interpolate data. Deterministic methods can either be exact or inexact interpolators. They don’t generate measures of estimate uncertainty. The stochastic interpolation methods used are Inverse Distance Weighting (IDW), Ordinary Kriging (OK), Universal Kriging (UK), and Delaunay Kriging (DK). Stochastic interpolation methods employed spatial correlations between values at neighboring points. The Kriging interpolation method used is regarded as the best linear unbiased estimator because its estimates are linear 35 combinations of the data with weights from non‐bias constraints. Kriging is computationally intensive when applied to large datasets, and it allows estimation of interpolation uncertainties. A semi‐variogram estimate was also used with isotropic models such as exponential, Gaussian, spherical, power law, linear, spartan and marten models. The results show that IDW presented more accurate for the basin than other values. DK provides the best cross‐ validation estimate for the extreme low value due to its localised nature of interpolation. Geostatistical technique in modeling uncertainties Geostatistics is an interpolation method used to estimate parameter values in space where there are no available data. It uses methods like kriging, co‐ kriging, variogram and trends to estimate parameters. Lin et al. (2017) applied geostatistical approach and Modflow to generate spatial distribution of hydraulic conductivities of a multiple aquifer in Pingtung, Taiwan. Geo‐statistic and stochastic elements work extremely together. This fact was discovered in its applications on many earth science studies. However, not all stochastic models are linked to geo‐statistics, leading to spatio‐stochastic modeling. Spatio‐stochastic models predict the behavior of complex systems with explicit knowledge of limited parameters relates to spatial relationships among locations in a system (Fortin et al., 2003). Geo‐statistics uses simulation techniques such as Lower Upper (LU) decomposition, turning bands, Sequential Gaussian (SG) and Simulated Annealing (SA) to model uncertainty of a regional variable e.g. hydraulic conductivity. Lin et al. (2017) used Simulated Annealing simulation to disturb pixel grid over the numerical iterations until pixel values reach the given histogram and semi‐variogram model. 36 They also combined Generalized Likelihood Uncertainty Estimation (GLUE), semi‐variogram model and simulation algorithm to generate random function model. GLUE is a kind of extended Monte Carlo random simulation that quantifies uncertainty. GLUE estimate likelihood of all possible outcomes for a specific distribution of inputs and determines behavioral parameter sets of a model. However, there have been critics of GLUE method for not being formally Bayesian. The weakness of geo‐statistical technique is that only a limited number of realizations can be generated making it difficult to characterize spatial uncertainty. Inverse method for unknown Boundary Conditions (BC) In case of unknown boundary conditions, a study carried by Irsa and Zhang (2012) proved that through the new inverse method, boundary conditions are not always needed for estimating aquifer parameters such as hydraulic conductivity. The Inverse method with a set of hybrid formulations was applied to the systems with problems of regular and irregular geometries, different heterogeneity patterns, variances and error magnitudes, and yielded successful results. BC are generally unknown in real aquifers and when wrong parameters are assigned, the parameter results estimated will be inaccurate and non‐unique leading to poor quality models. Such cases usually occurred when data used contain errors. Extended Kalman Filter (EKF) for parameter estimation in leaky aquifer Yeh and Huang (2005) used EKF method to identify hydraulic parameters in leaky aquifers systems both with and without considering aquifer storage. They found out that Kalman Filter is suitable for estimating parameters in a linear system; whilst Extended Kalman Filter is suitable for estimating hydraulic parameters in a nonlinear system. 37 Sequential Successive Linear Estimator (SSLE) SSLE technique was applied in a confined aquifer to solve transient and steady flow problems. SSLE inverse particularly solve transient flow problems, whilst Quasi Linear Bayesian Geostatic inverse method improve the distribution of depth averaged hydraulic conductivity. SSLE has been recognized as computational efficient, because the pumping or injection events are calculated separately and integrated by inner iteration loops, thus reducing the sizes of matrices for stochastic simulations (Yeh and Liu, 2000). 38 CHAPTER 4: DISTRIBUTION Distribution in general can be defined as a relationship for a sample space representing all possible outcomes of a random variable being observed in a particular system. The sample space meant here could be a set of real numbers, higher dimensional vector space or list of numerical or non‐ numerical values. Distributions are described in terms of their density or density functions. Density functions referred to how the proportion of data or likelihood of the proportion of observations changes over the range of the distribution. This kind of observation changes can be one of the reasons for applying distribution methods in estimating the hydraulic parameters stochastically. For density functions, the data should be classified as Probability Density Function (PDF), which calculates the probability of observing a given value, or Cumulative Density Function (CDF), which calculates the probability of an observation equal, or less than a value (Maus and Nijhawan, 2008). PDF and CDF are continuous, but the equivalent of a PDF for discrete distribution is termed as Probability Mass Function (PMF).Continuous distribution focuses on a set of ranges in a continuous way. Discrete is described by a set of values of a random variable. Walck (2007) stated that when analyzing raw data and fit the right distribution to the data, first figure out whether the distribution method and data can take discrete or continuous value, as some data may be restricted to discrete or continuous. In addition, it is advised to look at symmetry of the data and if there is asymmetry in which direction it lies if positive or negative outliers, equally likely or more likely than other. Analyze if there are upper or lower limits on the data because some data values can be less than zero while others 39 can exceed 100. The likelihood of observing extreme values in the distribution that may occur very infrequently or occur more often should be analyzed as well. It is often a challenge to choose the best distribution method suitable to estimate aquifer parameters, but in this investigation, we used mean, standard deviation, variance, tendency, skewness and kurtosis to estimate our aquifer parameters. It helps to figure out which distribution method better fits each specific parameter. Various distributions such as normal, log normal, Weibull, and others will be used to find the potential values for parameters, which are estimated to vary in space. 4.1 Statistical data analysis Statistical analysis describes the relationship between the data and their parameters by linking them in a statistical form. There are some hydraulic parameters e.g. storativity, leakage or transmissivity, which are sensitive, may not be known accurately. In this study, they are represented as randomly, sorted out by means of mean, variance, standard deviation, skewness, and kurtosis. The statistics of random variables will be solved by numerical equations. 4.1.1 Mean Mean can be arithmetic or geometric average mean of a dataset. Arithmetic mean can be defined as the sum of the values in a data set, divided by the count of number of data points. Calculating values using arithmetic mean is often regarded the best because it takes random variables and its independent, does not get affected by any poor performance of other values (Walck, 2007). The formula (15) for arithmetic mean: 40 1 (15) … = = Geometric mean can be calculated using equation (16) whereby you adds one to each value in a dataset to avoid problems that may arise from value with negative, and then multiply all the values together and raise the values to the power one (root square) by dividing by the count of number of values in a dataset. Geometric mean is best when calculating percentage changes in a datasets as it provide accurate measurements. The geometric mean formula: (16) = ⋯ 4.1.2 Variance Variance is calculated from the mean value whereby one has to first figure out the differences between the mean and average of the sample data, and square the values obtained from the sample differences to obtain the variance value. 4.1.3 Standard deviation It measure the amount of dispersion or variation of a sampled data. Here, one has to square root the value obtained in variance to obtain standard deviation. For the hydraulic parameters that may fit gaussian/normal distribution will be using standard deviation. The standard deviation formula (17) used to calculate sampled data: 41 (17) 1 = ( − ̅)² − 1 4.1.4 Skewness When data plotted on a graph is not centered (asymmetry), as like for normal distribution, but skewed to the left or right, resulting in one tail longer than the other, its skewness. Skewness can be calculated by using the average mean or standard deviation. A positive value of skewness indicates that mode is to the left of the mean and the longest tail is to the right of the bell shaped curve. Whilst a negative skewness shows that mode is to the right side of the mean, and longest tail of the bell shaped curve align to left side. When skewness is zero, the distribution is regarded as symmetry. 4.1.5 Kurtosis Kurtosis characterizes the tails of a distribution. The tail in a distribution with more extreme data values can be fatter indicating a positive kurtosis. When the tail is thinner due to less data values, the distribution will show a negative kurtosis. A distribution with zero kurtosis has the same outlier character as a normal distribution. 4.2Density distributions Probability density distribution is a function of a random variable taken from general sets of data, which falls within a particular range of a parameter values (Maus and Nijhawan, 2008). Any data that follows a probability density distribution is valuable and should be identified to which distribution it follows. There is no standard way of fitting data for distribution, we usually use 42 histograms to plot the data and then we analyze to which shape the histogram resembles, it can be Gaussian/normal, beta, or gamma distribution. This is because the rate of a function represent on the data in a form of physical quantities can change anytime. The aquifer parameters considered on this dissertation, all have distribution range values in aquifer; for example, transmissivity has a range value of 1 ‐ 100,000m2/d. For distribution, probability density function of a random variable can be either discrete or continuous, as discussed above. Continuous random variable is suitable for uniform distribution as it takes values in a continuous way. Non‐uniform distribution associates with discrete values. On this dissertation, we will analyze the data statistically to calculate the mean, variance, standard deviation, etc. to obtain average parameter values. The values obtained will be plotted on histogram and variogram graphs to determine which best probability distribution the data followed to estimate the aquifer parameters. Statistical software and manual plots can also help to plot the data on the graphs. The plots shows when the data follows a straight line on the graphs which indicates that the data follows the particular distribution, if it doesn’t fit, other distribution methods can used to assess the data. 4.2.1 Normal /Gaussian distribution Normal distribution, which is also known as Gaussian distribution is a kind of distribution which follows a symmetric shape. It can be obtained from the mean (µ) and standard deviation () or (s) of the distribution. So whenever the values of mean and standard deviation are changed, so does the positions and shapes of the distributions changes. 43 In a data set, we will identify which data follow a normal distribution or not. On a graph, normal distributions are characterized when plotting data, half of the values align to the left of the center and the other half to the right; the total area under the curve is 1; and the curve of the distribution is bell‐shaped and symmetrical. It’s easy to estimate the parameters on this kind of distribution (Maus and Nijhawan, 2008). The formula (18) used to calculate normal distribution: 1 () (18) () √2 Probabilities from different normal distribution can be converted to one probability involving the standard normal distribution in a process called standardization. 4.2.2 Binomial distribution It measures the probabilities of the number of success over a given number of chances tried with a specified probability of success in each try. Binomial distribution has only two possible outcomes: success and failure. The parameters of a binomial distribution are total number of samples and probability of success in each sample (Gleeson et al., 2016). 4.2.3 Lognormal distribution The lognormal distribution is the kind of distribution that shows asymmetrical and positive skewed shape. It does not take negative values in a dataset, only positive values. 44 4.2.4 Binomial–Exponential 2 distribution Bakouch et al. (2017) describes binomial exponential 2 as a distribution of a random sum of independent exponential random variables when sample size has a zero truncated binomial distribution. Stochastic models that accommodate zero value use this distribution type to observe different measures such as average, maximum and minimum in estimating distribution parameters. 4.2.5 Poisson distribution Gleeson et al. (2016) indicates that Poisson distribution measures the likelihood of a number of events occur in a given time interval to get the key parameter which is the average number of events in the given interval. It focuses on the number of occurrences of the event. Poisson distribution can be applied to estimate parameters in modeling homogenous type aquifers. Poisson distribution can represent the flow rate per time, which can vary. 4.2.6 Geometric distribution It focuses on the number of success in the number of trials by assuming that you were measuring the likelihood of when the first success will occur. The geometric distribution can be applied to water quality samples conducting the allowable parameters for human consumption from a certain water supplies. For aquifer pumping tests, geometric distribution can also be applied when conducting more than one‐step tests to determine the ideal discharge rate of a well (Freeze, 1975). 45 4.2.7 Exponential distribution Itis the probability distribution that describes the time between events.Exponential distribution focuses on the interval of time. During pumping test of a well, the time between when the pump was switched off, and the time when the groundwater water level recovered to its original static level (Freeze, 1975). 4.2.8 Generalized exponential distribution It shows an increasing or decreasing hazard rate depending on the shape parameter on a graph. It has many properties that are quite similar to those of the gamma distribution, but it has a distribution function like that of the Weibull distribution, which can be computed simply. The generalized exponential family has likelihood ratio ordering on the shape parameter; so it is possible to construct a uniformly most powerful test for testing a one‐sided hypothesis on the shape parameter when the scale and location parameters are known. 4.2.9Bernoulli distribution It is a discrete probability distribution that indicates a result from experiment or sampled data of what expected and actual outcome of the data results (Maus and Nijhawan, 2008). Just like binomial, there are also two possible outcomes for a Bernoulli distribution: success or failure. Binomial and Bernoulli have independent trails. An example of Bernoulli is investigating an aquifer to determine the allowable pumping rate, from a trial of random rates, of a well. 46 4.2.10 Weibull distribution The hydraulic parameters values in our work changes over time and Weibull distribution can handle this type of data. It is often preferred for analysing lifetime data in presence of censoring because it is much easier to handle than gamma distribution and it fits very well in many datasets as stated by Freeze (1975). We are using probability of Weibull distribution to help us to understand the groundwater flow behaviour at a time in space. The challenge is that the distribution of the sum of independent and identically distributed Weibull random variables is not simple to obtain. Therefore, the distribution of the mean of a random sample from a Weibull distribution is not easy to compute. 4.2.11 Beta distribution It is a continuous probability distribution having two parameters. It is applied mostly to random variables with limited ranges e.g. a – b, or 1-2. One of its most common uses is to model one's uncertainty about the probability of success of an experiment (Gleeson et al., 2016). 4.2.12 Gamma distribution It describesright‐skewed, continuous probability distributions. Its shape parameter is not confined to integer values. The gamma distribution starts at the origin and has a flexible shape. It is easy to estimate parameters by matching moments. 47 4.2.13 Rayleigh distribution Just like exponential distribution, Rayleigh also time dependent but it function linearly. Its reliability function decreases at a much higher rate than that of exponential distribution (Dey, Raheem and Mukherjee, 2017). 4.1.14 Cauchy distribution It has a probability distribution shape with much heavier tails than normal distribution. Its probability density is symmetric and unimodal. 4.1.15 Chi-squared distribution The chi‐squared distribution is obtained by adding up the squares of a number of normal variates (e.g. x2). Its distribution is used to test the goodness of fit between the observed data points and the predicted values by the model, which relate to the differences being normally distributed (Freeze, 1975). 4.1.16 Empirical distribution The empirical distribution is estimated directed from sample data without assuming the underlying algebraic form of the distribution model. 4.1.18 Pareto distribution The Perato distribution is a kind of distribution used to describe observed data. It is often model the tail of another distribution, either skewed to the other direction (Freeze, 1975). 4.1.19 Mittag-leffler distribution Mittag‐Leffler can take random values, which may be in a set of non‐negative or positive integers. It varies in shapes like a constant with unique mode at 48 zero, unimodal with one or two zero modes. Its distribution can be under, over or equi‐dispersed (Gleeson et al., 2016). 4.3 Estimation methods Below is the list of estimation methods used to estimate parameters in statistics (Bakouch et al., 2017): 1. Maximum Likelihood Estimation (MLE) is good in terms of consistency, efficiency, invariance and has good theoretical properties, but don’t perform well in small samples (Dey, Moala and Kumar, 2018) 2. Method of Moment estimation (MME) is obtained by comparing the first two theoretical moments with the sample moments. It is easily applicable and in some cases gives explicit forms for estimation of unknown parameters, but not for parameters of Weibull and Gompertz distributions. It also does not perform well in small samples. 3. Moments and pseudo-moments compares the first two theoretical moments of one with the sample moments. 4. Modified moments use the sample mean and variance. 5. Least-squares estimation (LSE) estimates the parameters of Beta distributions. Two variants of least squares methods can be used with the expectations and variances. 6. Percentile estimation (PE) can be used for cases like if the data come from a distribution function, which has a closed form, and the unknown parameters can be estimated by fitting straight line to the theoretical points obtained from the distribution function and the sample percentile 49 points. This method performs well in Weibull distribution and generalized exponential distribution. 7. L-Moments estimation (LME) compares the first two samples L‐ Moments with the corresponding population L‐Moments. This method has been used to estimate unknown parameters of generalized exponential distribution based on linear combination of order statistics. LME method is relatively strong to the effect of outliers whenever the mean of the distribution exists, as some higher moments might not exists making it more accurate in small samples. LME is efficient, compared to other estimation methods like MLE for normal distribution. 8. Maximum Product of spacing (MPS) and minimum spacing distance is an alternative to MLE for the estimation of parameters of continuous univariate distribution. 4.4 Hydraulic parameters The parameters of interest in this study are transmissivity, storativity and leakage factor. Many modelers treated these parameters as constant all time in a system, but it is not always valid due to variability of aquifer materials over time to time. We consider all these parameters being inconstant in aquifer system at a period. Transmissivity is obtained when hydraulic conductivity and thickness of a particular aquifer is known. Obtaining the hydraulic parameter values of an aquifer is usually sensitive because data used from observation wells at various distances taken at any location within the aquifer for example to estimate storativity, some people usually conclude it as the storativity of the whole domain of aquifer, which is inaccurate. 50 Not all aquifers are used for fresh water consumptions; some underground aquifers are used for disposal of wastewater and other hazardous liquids (Dagan and Neuman, 2005). Thus, it is very important to continuously monitor any leakage possibilities between aquifers for the best safety of groundwater resource. 4.4.1 Transmissivity Transmissivity (T) measures the rate of flow under a hydraulic gradient per unit through a cross section per unit width over the completely saturated thickness of a water‐bearing layer. It is expressed in m2/d. Transmissivity values are obtained from pumping test data with the use of Theis equation especially in confined aquifers, and Hantush Jacob in leaky aquifers. In situations when pumping test data are difficult and expensive to obtain, transmissivity can be derived from specific capacity data. Specific capacity data can be provided by drillers from step drawdown tests, which describe the performance of the well. The pumping rate, and discharge rate divided by drawdown is often used to get specific capacity data. Transmissivity affect the shape of drawdown on a cone of depression. If transmissivity is high, the drawdown will be broad and shallow, but when transmissivity is low, the cone of depression will indicate a narrow and deep drawdown shape. Estimation of transmissivity values is improved when using both pumping test (Theis) and specific capacity data (Theis, 1963). Transmissivity is calculated by using formula (19) integrating hydraulic conductivity over aquifer thickness. = (19) 51 In confined aquifer, for steady state flow Thiem equation (20) is used to calculate transmissivity: (20) = 2(ℎ − ℎ) For unsteady state flow, equation (21) can be used: (21) = () 4(ℎ − ℎ) For confined aquifers, distance drawdown, Cooper Jacob method (22) is used, as per one logarithmic cycle: 2.30 (22) = 4∆ Transmissivity depends on hydraulic conductivity and aquifer thickness making its values varies in different aquifers and in different places in the same aquifer. Among all hydraulic parameters, transmissivity is considered as higher than other parameters which varies in space, thus it should not be estimated as constant at all (van Leeuwen et al., 1998). Values estimated in principal aquifers ranges from less than 1 m2/d while in fractured aquifers ranges up to 100,000 m2/d for cavernous limestone and lava flows. In leaky aquifers, Hantush‐Jacob equation (23) calculates transmissivity: (23) = , 4 For storativity, equation (24) can be used: 52 4 (24) = ² where, Q: Pumping rate (m3/d) t: time since pumping began(days) r: distance from pumping well to observation well (m) B: leakage factor (m) 4.4.2 Storativity The capacity of how the water bearing materials are able to pass through and store water is termed as storativity. In confined aquifers, storativity occurs during compressibility of the water and the aquifer structure changes. During pumping, water comes from decompression of water and sediments. The aquifer materials are not drained as it remains saturated. Storativity is generally dimensionless, in confined aquifer it is the effect of specific storage and aquifer thickness. The equation (25) below used for calculating storativity in which units in the numerator and denominator cancel each other: (25) = = ( )( ℎ ℎ) ()() For confined aquifers, Cooper‐Jacob equation (26) for distance drawdown method calculates storativity as: 2.25 (26) = 53 The estimated value of storativity depends on whether the aquifer is confined, leaky or unconfined. Storativity values range from 1x10‐5 to 1x10‐3(0.001‐ 0.00001) for confined aquifers. The lower value of storativity indicates the deep drawdown expected. The higher value of storativity indicates that the drawdown will be minimal. 4.4.3 Leakage factor Leakage is a vertical transfer of water directed upward or downward through aquitard from overlying or underlying aquifers. Confined aquifers discharge water by leakage through its aquitard with a low hydraulic head. It’s very important to monitor and measure leakage to get the drawdown caused by pumping in an aquifer. Because some underground aquifers might be used for disposal of wastewater and other hazardous liquids or contamination from percolation of surface water. It is very important to continue monitoring any leakage possibility for the best safety of groundwater resource (Wang et al., 2018). Another reason for including leakage in groundwater monitoring is the permeability, which might change after earthquakes incident affecting the geologic fractures. Permeability in this case will likely to move not only in horizontal direction but also in vertical direction, as preexisting fractures are likely to be randomly oriented. The possibilities of leakage also exist at points where aquitard discontinued within subsurface. Leakage in an aquifer can be accessed with observation wells placed at reasonable location and depth to monitor drawdown. The changes in water level and recovery should also be monitored in a sufficient period. The drawdown data should be carefully corrected by means of tidal corrections or other methods applicable for analyzing the data to estimate leakage. When the 54 drawdown data has been corrected, it can analyzed with appropriate solution such as Hantush and Jacob (1955) solution, Boulton (1973) solution, Hunt and Scott (2007) solution, Zhan and Zlotnik (2002) solution, and Hunt stream (2003) and Hunt spring depletion (2004) solutions. Leakage factor determines the distribution of leakage through one or two aquitards into a leaky aquifer. The formula (27) is used to calculate leakage is: = √KHc (27) Formula (27) above can also be derived to the following (28): (28) m′mK = K The leakage has the dimension of length, usually expressed in meters as its influence is distributed over the area where leakage takes place. When analyzing the data, the large value of leakage factor shows a low leakage rate through the aquitard, whereas the small values indicates high leakage rate. High transmissivity values of the aquifer as well as high hydraulic resistance values of the aquitard leads to large values of leakage. Copty, Sarioglu and Findikakis (2006) assessed the impact of leakage on the transmissivity of the equivalent homogeneous leaky aquifer system using the value estimated for transmissivity as for arithmetic and geometric mean. As they relied on the assumption that head in un‐pumped aquifer cannot be affected by pumping as in pumped aquifer. It turned that the head in the un‐pumped aquifer unaffected by leakage results in larger estimate of transmissivity. 55 CHAPTER 5: ANALYSIS OF STOCHASTIC MODELS OF GROUNDWATER FLOW: CONFINED AND LEAKY AQUIFERS Stochastic approaches have been used in many fields of science, technology and engineering in the last decade with great success, as they are able to capture heterogeneities that cannot be described with normal deterministic approach. It is worth noting that several physical occurrences display statistical properties and this cannot be handled with normal deterministic equations as all the input parameters are assumed to be constant. In normal natural problems, one observed in normal situations that input parameters are stochastic not constant. This is the same situation with aquifer parameters including the transmissivity, storativity and hydraulic conductivity, they are geological dependent parameters and more importantly to the discharge rate which is the force applied by human. For instance, transmissivity cannot always statistically independent in the formation (space). The value of each specific parameter in space is independent as geological materials vary from one point to another depending on the fracture type. Fractures are made up of void space in rocks, which affect the flow properties and physical properties of a rock. Hydraulic parameters cannot be stationary in an entirely formation, like transmissivity can only be stationary at a sub unit of a geological formation in space, but not in the entirely formation. Boundaries, which might extend far away in some large‐scale aquifers, cannot be modeled well because of limited knowledge about them, so they should be treated as stochastic. If these variations are being neglected, it might lead to inaccurate results from the model produced misleading results to the decision makers to make recommendations for the future and it have negative impact 56 on the assessment of groundwater quantity e.g. average discharge or groundwater availability in an aquifer. It is very important to assess any uncertainty that may arise from the system. In this dissertation therefore, to be more representative, we shall consider the latest version of groundwater flow model within a confined and leaky aquifer proposed by Mathobo and Atangana (2018), and Amanda and Atangana (2018). The input parameters will be transformed into the stochastic one by assigning each parameter as a distribution. Without any doubt the transmissivity, leakage and storativity of the aquifer are not constant as we explained earlier, since one does not have a fixed function or law that these aquifer parameters follow, one will rely on stochastic approach to represent them. We shall start here with confined aquifer then leaky aquifer later. 5.1 Analysis of stochastic model of groundwater flow: Confined aquifers We analyzed the stochastic model for groundwater flow equation (29) in confined aquifer: (,) (,) ²(,) (29) = + ² where, T represents aquifer transmissivity, S for aquifer storativity in a (r) radial distance at a (t) time. T is considered as constant in Theis model. However, in reality T changes as a function of s (drawdown) as we argued above, in this case: ∈ (,,,,⋯,) 57 Drawdown data can be analyzed by using statistical methods such as arithmetic mean or geometric mean to obtain an accurate value represented by distribution probability techniques assigned in a stochastic model. With arithmetic mean, use equation (30): 1 (30) = Therefore geometric mean, we have the following (31): (31) = = ⋯ We consider0 ˂ ≤ 1, we call it stochastic parameter. We consider as a distribution associated to . can be normal, log normal, Pareto, Poisson, Mittag‐Leffler distribution, etc. We replace by T to obtain (32) = + (32) Following the same approach, we obtained (33): ̅ = ̅ + (33) Now our stochastic groundwater flow model in confined aquifer is given as equation (34): 58 (,) (,) ²(,) (34) = + ² The rate of drawdown changes in time per radial distance due to velocity or dispersion of water when distributed in a confined aquifer. This condition occurred during pumping at different discharge rates, the confining layers will moves down causing the aquifer to compress. The pumped water comes from storage within the aquifer at the beginning of pumping and drawdown increases resulting in unsteady state flow condition until such a time a well indicate a small change in drawdown (steady state) as storativity increases with the aquifer’s thickness. 5.2 Analysis of stochastic model of groundwater flow: Leaky aquifers Analyzing hydraulic parameters in leaky aquifers, it’s a matter of radial distance of a pumped well and observation well which needs to be placed farthest away apart to measure significant leakage in aquitard. Theoretically, at the beginning, the well is being supplied by the aquifer storage but at a later stage, the well will be supplied by storage from aquitard. The effect of leakage at this point will be minimal but needs to be considered as it also contributes to the drawdown. The value for leakage factor should be determined from drawdown of the pumped well and observation wells. Transmissivity and storativity values should be obtained using the drawdown data from pumped well only because effect of leakage is not significant and not influencing drawdown. The classical equation (35) commonly used in the literature is given below as: 59 ℎ ℎ ∆ 1ℎ ℎ(,) (35) = 1 + + + ƛ² Following the discussion presented for the case of confined aquifer, the above equation can be converted to the following equation (36) but with a new parameter introduced representing leakage factor: ℎ ℎ ∆ 1ℎ ℎ(,) (36) = 1 + + + ƛ² Both stochastic models cannot be solved using classical analytical methods as they become nonlinear due to stochastic coefficients. Therefore, where the analytical method failed, numerical methods can be used to provide an approximated solution. This will be presented in detail in the next chapter. 60 CHAPTER 6: APPLICATION OF NEWTON METHOD ON STOCHASTIC GROUNDWATER FLOW MODELS FOR CONFINED AND LEAKY AQUIFERS The Newton method has been regarded as the very efficient in solving nonlinear equations (Burden and Faires, 2005) such as Partial Differential Equations (PDE) mostly used in this dissertation. Differential equations, which can be partial or ordinary, describe the relationship between the physical quantities and their rate of change, which in this case known as derivative. The rule of Newton method for initial guess to solve equation is ()= 0. In general, it is used to approximate and solve equations numerically by using iteration method, discretization, Taylor series or Jacobian method. Solving equations by using iterative process, you have to repeat the procedure repeatedly to find the root of an equation. Taylor series is achieved by performing a series of expansion on a function of that point. In this chapter, we are using Finite Difference to solve PDE using discretization to reduce data size, which corresponds with temporal (time) and spatial (space) and applying the Newton’s method to approximate an exact solution of the parameters, which will be assigned in a model of Theis stochastic for confined aquifers and Hantush stochastic for leaky aquifers. Transmissivity is measured at a time, and storativity is measured per radial distance in space, which all varies. Unlike Ordinary Difference Equation (ODE) which only requires certain ordinary information (one variable), for PDE which make up most of the equation used on this dissertation, we need to know initial values which may contain multiple variables and extra information about the behavior of the solution at temporal and spatial domain. We noticed that the best way to achieve all that we can use the process of discretization to get the exact approximate for the solution. 61 We begin with the basic of derivatives with their functions, mainly the 1storder derivative of a function, the 2nd order derivative of a function, and Crank‐ Nicholson to prove their accuracy so that we can apply them into our model. We also include the forward, backward and central difference techniques to approximate the derivative of functions. For 1storder derivative with functions with one or two variables is given as equation (37) and (38): ( + ℎ)− () (37) ′()= lim → ℎ ( + ℎ)− () (38) ()= ()= lim → ℎ For 2ndorder derivative with functions of one or two variables is given as equation (39) and (40): () ( + ℎ)− () (39) ()= ′() = = lim → ℎ ()() ()() lim − = → ℎ (40) ( + ℎ)− 2()+ ( − ℎ) ′′()= lim → ℎ² Forward and backward difference for approximating derivative of functions with one and multiple variables indicated on equation (41) to (43): ( )− () (41) ()= ∆ 62 ()− ( ) ()= ∆ (,) (,)− (,) (42) = ∆ (,) (,)− (,) (43) = ∆ In case where there is a recharge to the aquifer, a sum of Forward and Backward difference is obtained to get what is called Central difference in equation (44): ()− () (44) ′()= ∆ where, ∆: step site : point of interest : point of interest plus one : point of interest minus one Approximate of a solution with2ndorder derivative for a function with one and two variables using forward, backward and central difference method from equation (45) to (46): 63 ()− 2()+ () (45) ′′(ᵢ)= ∆² ²(,) (,)− 2(,)+ (,) (46) = ² (∆)² Crank Nicholson (equation 47): ²(,) 1 (47) = [(,)− 2(,)+ (,)² (∆)² + (,)− 2(,)+ (,)] Introducing the Newton method, we find the following starting from equation (48 to (52) : () (48) = [,()] ()− () 1 (49) = [(,())+ ( , )] ∆ 2 − 1 (50) = [( , )+ ( , )] ∆ 2 ∆ (51) = + [(, + ( , )] 2 1 − (52) − ∆[(, )+ (, )] = 0 2 64 From equation described above, now we replace = , and = () equation (53) to (58): ∆ (53) − ()− [(, )+ (,)] = 0 2 ∆ (54) ()= − − [(, )+ (,)] = 0 2 () ∆ (,) (55) ′()= = 1 − 2 ∆ (,)¹ (56) ′()= 1 − 2 () (57) = − () ∆ − [(, )+ ( , )] (58) = − ∆ ( 1 − , ) In the next section, we approximate the solution for stochastic confined and leaky aquifer models using PDE by means of discretization of time and space. We first start with confined aquifer, and thereafter leaky aquifer as derived in the previous chapter. 65 6.1 Application of Newton method on stochastic Theis’s confined aquifer To apply Newton method, we reformulate the groundwater flow equation in confined aquifer to achieve a stochastic model, which covers the rate of change in transmissivity at a time, and storativity per radial distance as follow in equation (59): (,) ²(,) (59) + = (,,(,)) ² At the points of (,), werearrange to obtain equation (60) until (61) : (,)− (,) (60) ∆ 1 = [(,,(,))+ (,),( , )] 2 ∆ (61) − = [(,, )+ (,, )] 2 From the previous equation, we then set = = to get the root of the function and it gives us the following equation (62), (63): () ∆ (62) = + [(,, )+ (,,)] 2 ∆ (63) − − [(,, 2 )+ (,,)] = 0 66 According to the Newton method, it should be expressed as equation (64) to (66): ∆ (64) ()= − − [( , , )+ (,2 ,)] () ∆ (,,) (65) = 1 − 2 ∆ (,,) = 1 − 2 ∆ ∆ − , , − ( , , ) (66) − = ∆ ( , , ) 1 − To further discretize the domain, we recall that equation (67) to (70): (,) (,) (67) (,,(,))= + 1 − − 2 + (68) ,,(,) = + ∆ (∆)² 1 − ( ) − 2 + (69) ,,(,) = + ∆ (∆) ( ,,) 2 2 (70) = − = − (∆) (∆) 67 And the final product after all the steps above, our groundwater flow stochastic equation for confined aquifer with the application of new scheme provided the following solution, equation (71). The solution proves that the rate of drawdown changes per radial distance in respect of dispersion of water at a time in an aquifer. The discharge rate at which the well is being pumped also affects the drawdown. = ∆ + + + ∆ (∆)² ∆ (∆) − 2 ∆ 1 − (∆) (71) 6.2: Application of Newton method on stochastic Hantush’s leaky aquifer In order to achieve the groundwater flow stochastic model for leaky aquifer, we also applying the Newton method just like the similar procedure we did for confined aquifer. We therefore reformulate groundwater flow equation of leaky aquifer, and apply the Newton method to approximate the solution (72): ²(,) (,) (72) + + = (,,(,)) ² ƛ² At the points of (,), we obtain (73) and (74): 68 − 1 (73) = [(,, )+ (, )] ∆ 2 ∆ (74) − = [( ,, )+ (, , )] 2 Since = and if = , then the solution will be (75) and (76): () ∆ (75) = + [(, , )+ (,2 ,)] ∆ (76) − − [(,, )+ (,,)] = 0 2 By generalizing the above equation, we apply the Newton method helps us to obtain the following (77), (78) and (79): ∆ (77) ( )= − − [(,, )+ (,,)] 2 This implies that: (78) () ∆ (,,) = 1 − 2 ∆ (,,) = 1 − − 2 ƛ² 69 ∆ ∆ − ,, − ( , , (79) ) − = ∆ (,,)1 − We remember the previous equation, and we get solution (80), (81), (82) and (83): (,) (,) (80) (,,(,))= + 1 ( ) − (,) − 2 + (81) ,, = + + ∆ ƛ² (∆)² ,,( ) (82) 1 − (,) = + ∆ ƛ − 2 + + (∆)² (,, ) 2 2 (83) = − = − (∆) (∆) After we discretized the stochastic groundwater flow for leaky equation with the steps specified above, unlike for confined aquifer’s equation, here we now 70 adding a new parameter ƛ representing leakage factor in an aquifer in equation (84), where: (84) = + ∆ + + + ∆ (∆)² ∆ (∆) ƛ − 2 ∆ 1 − + (∆) ƛ The new parameter is known as leakage factor, which shall be now considered even though its effects on drawdown are negligible. Therefore, the modified new scheme is made up of certain parameters in a leaky aquifer such as transmissivity, storativity, leakage factor and drawdown. Leakage factor has been ignored in Hantush – Jacob’s equation, thus we then modified it and include the leakage factor, which contribute to the storativity of the aquifer irrespective of the thickness of the aquitard. The fact that leakage factor should be recognized further balance the steady condition when the water leaking through the aquitard into the aquifer ideal for recharging the aquifer. 6.3 Stability for the stochastic confined aquifer equation Stabilizing the model means that most of the objects used in the model will be stable over time and will not need changes. The model is considered accurate when it reaches stability convergence and there is little difference between exact and approximate of a solution when fitting the dataset. If the model’s parameters are different during the forecast period than they were during the construction period, then the model estimated will not be very useful regardless of how well it was estimated (Sinha and Prasad, 1979). If the 71 model’s parameters were unstable over the construction period, then the model was not even a good representation of how the parameters actually occur within groundwater flow in the aquifer at a time. We analyzed our model for any inconsistency by means of stability before the predictions can be applied on a real system. We are now doing stability confined aquifer as stipulated in equation (85), (86) (87) to obtain (88) and (89): 2 1 + ∆ (85) (∆) 1 − (∆) ∆ = 1 − ∆ (∆) 1 − (∆) − ∆ 1 1 + + 2 ∆ (∆) ∆ 1 + (∆) ∆ 1 1 + + 2∆ (∆) ∆1 + (∆) − ∆ ∆ 1 + − 2 2(∆) ∆ 1 + (∆) − ∆ ∆ 1 + − 2 2(∆) ∆ 1 + (∆) If = ,then (∆) = 72 2 1 + (86) ∆ (∆) 1 − (∆) ∆ = 1 − ∆(∆) 1 + (∆) − ∆ 1 1 + (∆)⎛ + ⎞2∆ (∆) ∆ 1 + ⎝ (∆)⎠ ∆ 1 1 + (∆)⎛ + ⎞ 2∆ (∆) ∆1 + ⎝ (∆) ⎠ −∆ ∆ 1 + (∆)⎛ − ⎞2 2(∆) ∆ 1 + ⎝ (∆) ⎠ −∆ ∆ 1 + (∆)⎛ ⎞ − 2 2(∆) ∆ 1 + ⎝ (∆)⎠ We further convert: (∆) = . ∆ (∆) = . ∆ And add it to the solution to get: 73 2 1 + (87) ∆ (∆) 1 − (∆) ∆ = 1 − ∆ (∆) 1 − (∆) − ∆ 1 1 + . ∆ ⎛ + ⎞2∆ (∆) ∆ 1 + ⎝ (∆) ⎠ ∆ 1 1 + . ∆ ⎛ + ⎞2∆ (∆) ∆1 + ⎝ (∆) ⎠ −∆ ∆ 1 + . ∆ ⎛ − ⎞2 2(∆) ∆ 1 + (∆)⎝ ⎠ −∆ ∆ 1 + . ∆ ⎛ − ⎞ 2 2(∆) ∆ 1 + ⎝ (∆)⎠ After simplification, we obtained: 74 (88) 2 1 + ∆ (∆) 1 − (∆) ∆ = 1 − ∆ (∆) 1 − (∆) − ∆ 1 1 + ∆ ⎛ + ⎞2 ∆ (∆) ∆ 1 + ⎝ (∆)⎠ ∆ 1 1 + ∆ ⎛ + ⎞2∆ (∆) ∆1 + ⎝ (∆) ⎠ −∆ ∆ 1 + ∆ ⎛ − ⎞2 2(∆) ∆ 1 + (∆)⎝ ⎠ −∆ ∆ 1 + ∆ ⎛ − ⎞ 2 2(∆) ∆ 1 + ⎝ (∆) ⎠ 75 2 ∆ 1 + − ∆ 1 − (89) ∆ ∆ (∆) 1 − (∆) 1 + (∆) (∆) − ∆ 1 1 + ∆ + 2∆ (∆) ∆ 1 + (∆) ∆ 1 1 = + 2∆ (∆) ∆1 + (∆) − ∆ ∆ 1 + ∆ ⎛ − ⎞ 2 ∆ 2(∆) 1 + (∆) ⎝ ⎠ −∆ ∆ 1 + ∆ ⎛ − ⎞ 2 2(∆) ∆ 1 + ⎝ (∆) ⎠ For simplicity, we put the following in equation (90) to get the outputs in equation (91), (92) and (93): 76 2 (90) = 1 + ∆ (∆) 1 − (∆) ∆ = 1 − ∆ (∆) 1 + (∆) − ∆ 1 1 = + 2∆ (∆) ∆ 1 + (∆) ∆ 1 1 = + 2∆ (∆) ∆1 + (∆) − ∆ ∆ 1 = − 2 2(∆) ∆ 1 + (∆) 2 = 1 + ∆ (∆) 1 − (∆) ∆ = 1 − ∆ (∆) 1 + (∆) − ∆ 1 1 = + 2∆ (∆) ∆ 1 + (∆) ∆ 1 1 = + 2∆ (∆) ∆1 + (∆) − ∆ ∆ 1 = − 2 2(∆) ∆ 1 + (∆) Our equation now becomes as: 77 a ∆ − + ∆ = a + ∆ + ∆ (91) [ + (cos∆ (92) + sin(∆))+ (cos ∆ − sin ∆)] = [ + 2 cos(∆)] [ + cos(∆) (93) + cos(∆) + sin(∆)( − )= ( + 2 cos(∆))] Factorizing the previous equation (93) to get equation (94): [ + ( + )cos(∆)+ sin(∆)( − )] (94) = [( + 2 cos(∆))] Every model is regarded stable when it converges this way: < 1 The approximate and exact solution for the real and imaginary system proved that our model is stable as it resulted in a solution gives< 1, which is accurate as theoretical predicted as proved in equation (95), (96), (97) and (98): + 2 cos(∆) (95) = + ( + )cos(∆)+ sin(∆)( − ) + 2 cos(∆) (96) = + ( + )cos(∆)+ sin(∆)( − ) 78 | + 2 cos(∆)| (97) = ( + ( + )cos(∆)) + (sin(∆)( − )) < 1 | + 2 cos(∆)| (98) < ( + ( + )cos(∆)) + (sin(∆)( − )) The numerical analysis for our stochastic equation in confined aquifer converges with the actual analysis of a system. The parameters assigned to construct the models meets with what theoretical predicted. When we arrived at this stage, it shows that our model will perform well and can be relied on to make prediction in the future. Therefore, in the next part of simulation, we should expect these parameters to be similar to the simulated ones. 6.4 Stability for the stochastic leaky aquifer equation Referring to the leaky equation, we follow the similar approach, as done for the confined aquifer, for stability. We should take note that for the leaky aquifer equation, there is an additional parameter known as leakage factor, which compromise the leaky’s aquitard, but not common in confined aquifer. This can be checked in equation (99) (100) (101) to obtain (102) and (103): 79 2 1 + ∆ (99) (∆) 1 − + (∆) ƛ ∆ = 1 − ∆ (∆) 1 − + (∆) ƛ − ∆ 1 1 + + 2∆ (∆) ∆ 1 + + (∆) ƛ ∆ 1 1 + + 2∆ (∆) ∆ 1 + + (∆) ƛ − ∆ ∆ 1 + − 2 2(∆) ∆ 1 + + (∆) ƛ − ∆ ∆ 1 + − 2 2(∆) ∆ 1 + + (∆) ƛ If = ,then = (∆) 80 2 1 + (100) ∆ (∆) 1 − + (∆) ƛ ∆ = 1 − ∆ (∆) 1 + + (∆) ƛ − ∆ 1 1 + (∆)⎛ ⎞ + 2 ∆ (∆) ∆ 1 + + ⎝ (∆) ƛ⎠ ∆ 1 1 + (∆)⎛ ⎞ + 2∆ (∆) ∆ 1 + + (∆) ƛ⎝ ⎠ −∆ ∆ 1 + (∆)⎛ − ⎞ 2 ∆ 2(∆) 1 + + ⎝ (∆) ƛ⎠ −∆ ∆ 1 + (∆)⎛ − ⎞ 2 2(∆) ∆ 1 + + ⎝ (∆) ƛ ⎠ We further convert: (∆) = . ∆ (∆) = . ∆ Moreover, add it to the solution to get: 81 1 2 1 + ƛ (101) ∆ (∆) 1 − + (∆) ƛ ∆ = 1 − ∆ (∆) 1 − + (∆) ƛ − ∆ 1 1 + . ∆ ⎛ + ⎞ 2 ∆ (∆) ∆ 1 + + ⎝ (∆) ƛ⎠ ∆ 1 1 + . ∆ ⎛ + ⎞2∆ (∆) ∆ 1 + + ⎝ (∆) ƛ⎠ −∆ ∆ 1 + . ∆ ⎛ ⎞ − 2 2(∆) ∆ 1 + + ⎝ (∆) ƛ⎠ −∆ ∆ 1 + . ∆ ⎛ − ⎞ 2 2(∆) ∆ 1 + + ⎝ (∆) ƛ⎠ After simplification, we obtained the following: 82 1 2 1 + ƛ ∆ (∆) 1 − + (∆) ƛ (102) ∆ = 1 − ∆ (∆) 1 − + (∆) ƛ − ∆ 1 1 + ∆ ⎛ + ⎞ 2∆ (∆) ∆ 1 + + ⎝ (∆) ƛ⎠ ∆ 1 1 + ∆ ⎛ + ⎞2∆ (∆) ∆ 1 + + ⎝ (∆) ƛ ⎠ −∆ ∆ 1 + ∆ ⎛ − ⎞2 2(∆) ∆ 1 + + (∆) ƛ⎝ ⎠ −∆ ∆ 1 + ∆ ⎛ − ⎞ 2 2(∆) ∆ 1 + + ⎝ (∆) ƛ⎠ 83 1 2 1 + (103) ƛ ∆ (∆) 1 − + (∆) ƛ ∆ − ∆ 1 − ∆ (∆) 1 + + (∆) ƛ − ∆ 1 1 + ∆ + 2 ∆ (∆) ∆ 1 + + (∆) ƛ ∆ 1 1 = + 2∆ (∆) ∆ 1 + + (∆) ƛ − ∆ ∆ 1 + ∆ ⎛ − ⎞ 2 2(∆) ∆ 1 + + (∆) ƛ ⎝ ⎠ −∆ ∆ 1 + ∆ ⎛ − ⎞ 2 2(∆) ∆ 1 + + ⎝ (∆) ƛ ⎠ For simplicity, we put the following parameters to obtain equation (104): 1 2 (104) = 1 + ƛ (∆) ∆ 1 − + (∆) ƛ ∆ = 1 − ∆ (∆) 1 + + (∆) ƛ − ∆ 1 1 = + 2 ∆ (∆) ∆ 1 + + (∆) ƛ ∆ 1 1 = + 2∆ (∆) ∆ 1 + + (∆) ƛ 84 −∆ ∆ 1 = − 2 2(∆) ∆ 1 + + (∆) ƛ Our equation now becomes as equation (105), (06) and (107): a − ∆ + ∆ = a + ∆ + ∆ (105) [ + (cos∆ (106) + sin(∆))+ (cos ∆ − sin ∆)] = [ + 2 cos(∆)] [ + cos(∆) (107) + cos(∆) + sin(∆)( − )= ( + 2 cos(∆))] Factorizing the previous equation (107) to get (108): (108) [ + ( + )cos(∆)+ sin(∆)( − )] = [( + 2 cos(∆))] The model is regarded stable when it converges this way, following the steps from equation (109 until last step (112) : < 1 + 2 cos(∆) (109) = + ( + )cos(∆)+ sin(∆)( − ) 85 + 2 cos(∆) (110) = + ( + )cos(∆)+ sin(∆)( − ) | + 2 cos(∆)| (111) = ( + ( + )cos(∆)) + (sin(∆)( − )) < 1 | + 2 cos(∆)| (112) < ( + ( + )cos(∆)) + (sin(∆)( − )) The numerical analysis for our stochastic equation in leaky aquifer converges with the actual analysis of a system. It means that the model parameters assigned during construction meets exactly with the theoretical predicted. On the next section of simulation, we should expect the model to perform well so that it can be used for predictions in the future. 6.5 Simulation After we developed our stochastic groundwater flow models in confined and leaky aquifers, the next step was to run the model via simulation to obtain the actual representative of the behavior of the system. Simulation helps to test the system and its condition being investigated which may be difficult to imitate in the real world. For complex systems, simulation provide valuable solutions by giving clear insights which is necessary for prediction or forecasting the future behavior of a system and also determine what can be done to influence that future behavior. Below are the images obtained after 86 the completion of simulation representing the groundwater flow behavior in confined and leaky aquifer. We consider a theoretical sample of storativity and transmissivity, we compute harmonic, geometric and arithmetic mean and we consider the normal distribution in this case. Numerical simulation are depicted in figure 2, 3 ,4, 5, 6,7, 8,9 and 10 including contour plots for confined case. The drawdown is plotted for this case. Figure 6.1 Numerical solution with harmonic mean for transmissivity and storativity 87 Figure 6.2: Contour plot for harmonic mean 88 Figure 6.3: Numerical simulation with arithmetic mean Figure 6.4: Contour plot with arithmetic mean 89 Figure 6.5: Numerical solution with geometric Figure 6.6: Contour plot of geometric mean 90 Figure 6.7: Numerical solution with geometric mean Figure 6.8: Contour plot with geometric mean 91 Figure 6.9: Numerical simulation with arithmetic mean Figure 6.10: Contour plot with arithmetic mean 92 Here we present the reduction of water level for the leaky case, numerical simulation are depicted in figure 11, 12, 13, 14, 15, 16, 17, 18. Figure 6.11: Numerical simulation for harmonic mean Figure 6.12: Contour plot with harmonic mean 93 Figure 6.13: Numerical simulation with harmonic mean Figure 6.14: Contour plot with harmonic mean 94 Figure 6.15: Numerical simulation for arithmetic mean Figure 6.16: Contour plot for arithmetic mean 95 Figure 6.17: Numerical simulation for geometric mean Figure 6.18: Contour plot with geometric mean 96 The numerical simulation show that, modelling with stochastic approach is more suitable as one may have many possibilities to capture heterogeneities. One of the most important finding in this work is the fact that the well‐known arithmetic mean that is widely used in the field of Geohydrology provides exaggerated information. This exaggeration can be well explained with the concept of salaries, the richest man leaving in same room with a poor man, if one asks what the average salary in this room is obviously one will conclude that each person in the room is rich if the average is arithmetic. Another illustrative example is that of quantity of water within an aquifer, one will drill four bore wells and get the quantity of water in each and then use the arithmetic average to say the average water within the portion is X for instance, while there can be some place within that system with no water. On the other hand, one will see that modelling using the harmonic mean will lead to underestimated results; however, one could get moderate results using the geometric mean, which is comprised between arithmetic mean and harmonic mean. While we suggest the use of stochastic models in Geohydrology to capture more heterogeneity, we should point out that, deterministic could be useful if the situation is purely homogeneous, of course this situation is very hard to be found in real life. 97 CONCLUSION The aquifers are being described by their lithology, geometry, flow behavior, and fractures. Confined and leaky are the aquifers considered on this dissertation. It is important to recall that, the flow of subsurface water is affected by the geometry of the aquifer within which the flow takes place. Therefore, assess the effect these aquifers to the groundwater flow, we estimated the aquifer parameters of which we assumed vary in space at a time. The aquifer parameters help to understand the behavior groundwater flow within the subsurface, and this can only be achieved by means of modeling. The modeling techniques such as stochastic and deterministic are commonly used to estimate aquifer parameters in groundwater flow system. Deterministic techniques can be used for small scale projects like in homogenous environment because they predict a single outcome. Stochastic technique can handle complicated projects in both homogenous and heterogeneous environment as ignored earlier. The latest version (Atangana and Ramatsho) is then modified to a new scheme indicating that aquifer parameters changes in time at a point in a system. Based on the latest version equation of Atangana and Ramatsho, and also Atangana and Mathobo, the mathematical equation was modified further to capture the statistical analysis setting of aquifer parameters such transmissivity, storativity and leakage factor. The statistical analysis including the mean, variance, skewness and standard deviation were used in the modified equation to estimate aquifer parameters. The numerical analysis of stochastic models for the confined and leaky equation was solved with application of the Newton method by means of discretization. Different methods of probability distribution were applied on the models to describe 98 appropriate better fit for each parameter. The model results show a good balance which means that stochastic technique captured well prediction of the system. When comparing with other simulation methods such as deterministic models done in groundwater system their strength relies on capturing any uncertainty, arise in any environment. Stochastic models produced in this dissertation contain PDE equations with functions of one and multiple variables which are analyzed numerically. The existing models of Theis, and Hantush and Jacob were reviewed and modified to a new version (Atangana and Ramatsho) that includes the geological matrix, which was, stochastic models for confined and leaky aquifers obtained for this dissertation are unique, accurate and good quality, ready to be used for current and future predictions. The numerical simulation suggests the use of stochastic while the geometric mean is perform. 99 REFERENCES Adem, G. and Batelaan, O., 2006. Modeling groundwater‐surface water interaction by coupling MODFLOW with WetSpa. Geophysical Research, 8, p.03181. Adhikary, P.P. and Dash, C.J., 2017. Comparison of deterministic and stochastic methods to predict spatial variation of groundwater depth. Applied Water Science, 7(1), pp.339‐348. Amanda, R. and Atangana, A., 2018. Derivation of a groundwater flow model within leaky and self‐similar aquifers: Beyond Hantush model. Chaos, Solitons & Fractals, 116, pp.414‐423. Bakouch, H.S., Dey, S., Ramos, P.L. and Louzada, F., 2017. Binomial‐exponential 2 distribution: Different estimation methods with weather applications. TEMA (São Carlos), 18(2), pp.233‐251 Burden, R.L. and Faires, J.D., 2005. Numerical Analysis 8th Edition. Thompson Brooks/Cole. Copty, N.K., Sarioglu, M.S. and Findikakis, A.N., 2006. Equivalent transmissivity of heterogeneous leaky aquifers for steady state radial flow. Water resources research, 42(4). Dagan, G. and Neuman, S.P. eds., 2005. Subsurface flow and transport: a stochastic approach. Cambridge University Press. Dagan, G., 1989. Flow and transport in porous media Dey, S., Moala, F.A. and Kumar, D., 2018. Statistical properties and different methods of estimation of Gompertz distribution with application. Journal of Statistics and Management Systems, 21(5), pp.839‐876. 100 Dey, S., Raheem, E. and Mukherjee, S., 2017. Statistical properties and different methods of estimation of transmuted Rayleigh distribution. Revista Colombiana de Estadística, 40(1), pp.165‐203. Fetter, C.W., 2001. Properties of aquifers. Applied hydrogeology, pp.66‐112. Fortin, M.J., Boots, B., Csillag, F. and Remmel, T.K., 2003. On the role of spatial stochastic models in understanding landscape indices in ecology. Oikos, 102(1), pp.203‐212. Freeze, R.A., 1975. A stochastic‐conceptual analysis of one‐dimensional groundwater flow in nonuniform homogeneous media. Water Resources Research, 11(5), pp.725‐741. Gelhar, L.W., 1993. Stochastic subsurface hydrology. Prentice‐Hall. Gleeson, T., Befus, K.M., Jasechko, S., Luijendijk, E. and Cardenas, M.B., 2016. The global volume and distribution of modern groundwater. Nature Geoscience, 9(2), p.161. Hantush, M.S. and Jacob, C.E., 1955. Non‐steady radial flow in an infinite leaky aquifer. Eos, Transactions American Geophysical Union, 36(1), pp.95‐100. Hantush, M.S., 1960. Modification of the theory of leaky aquifers. Journal of Geophysical Research, 65(11), pp.3713‐3725. Hsu, K.C., Zhang, D. and Neuman, S.P., 1996. Higher‐order effects on flow and transport in randomly heterogeneous porous media. Water Resources Research, 32(3), pp.571‐582. Hu, L., Chen, C. and Chen, X., 2011. Simulation of groundwater flow within observation boreholes for confined aquifers. Journal of hydrology, 398(1‐2), pp.101‐108. Irsa, J. and Zhang, Y., 2012. A direct method of parameter estimation for steady state flow in heterogeneous aquifers with unknown boundary conditions. Water Resources Research, 48(9). 101 Koch, M. and Arlai, P., 2007. Deterministic and stochastic Modeling of Groundwater Flow and Solute Transport in the heavily‐stressed Bangkok coastal Aquifer. Thailand, and investigation of optimal management strategies for possible aquifer restoration, IAH, pp.17‐21. Kruseman, G.P., De Ridder, N.A. and Verweij, J.M., 1970. Analysis and evaluation of pumping test data (Vol. 11). The Netherlands: International institute for land reclamation and improvement. Kumba, F.F., 2003. Farmer participation in agricultural research and extension service in Namibia. Journal of International Agricultural and Extension Education, 10(3), pp.47‐55. Lavallée, D. and Archuleta, R.J., 2003. Stochastic modeling of slip spatial complexities for the 1979 Imperial Valley, California, earthquake. Geophysical Research Letters, 30(5). Lin, Y.P., Chen, Y.W., Chang, L.C., Yeh, M.S., Huang, G.H. and Petway, J., 2017. Groundwater Simulations and Uncertainty Analysis Using MODFLOW and Geostatistical Approach with Conditioning Multi‐Aquifer Spatial Covariance. Water, 9(3), p.164. Mathobo, M. and Atangana, A., 2018. Analysis of exact groundwater model within a confined aquifer: New proposed model beyond the Theis equation. The European Physical Journal Plus, 133(10), p.415. Maus, G. and Nijhawan, R., 2008. Motion into and out of the blind spot: Evidence for spatial extrapolation of moving objects. Perception, 37(2), pp.310‐ 310. Neuman, S.P. and Witherspoon, P.A., 1969. Theory of flow in a confined two aquifer system. Water Resources Research, 5(4), pp.803‐816. Payne, S.M. and Woessner, W.W., 2010. An Aquifer Classification System and Geographical Information System‐Based Analysis Tool for Watershed 102 Managers in the Western US 1. JAWRA Journal of the American Water Resources Association, 46(5), pp.1003‐1023. Renard, P., 2007. Stochastic hydrogeology: what professionals really need?Groundwater, 45(5), pp.531‐541. ŞEN, Z., 2000. Non‐Darcian groundwater flow in leaky aquifers. Hydrological sciences journal, 45(4), pp.595‐606. Shvidler, M.I., 1964. Filtration flows in heterogeneous media. Consultants Bureau. Sinha, N.K. and Prasad, T., 1979. Some stochastic modelling techniques and their applications. Applied Mathematical Modelling, 3(1), pp.2‐6. Theis, C.V., 1963. Estimating the transmissibility of aquifer from the specific capacity of wells. Methods of determining permeability, transmissibility, and drawdown. US Geological Survey Water Supply Pater, 1536, pp.332‐336. Thiem, G., 1906. Hydrologic methods. Gebhardt, JM Leipzig, Germany. van Leeuwen, M., te Stroet, C.B., Butler, A.P. and Tompkins, J.A., 1998. Stochastic determination of well capture zones. Water Resources Research, 34(9), pp.2215‐2223. Varouchakis, Ε.A. and Hristopulos, D.T., 2013. Comparison of stochastic and deterministic methods for mapping groundwater level spatial variability in sparsely monitored basins. Environmental monitoring and assessment, 185(1), pp.1‐19. Wang, Y., Ocampo‐Martinez, C. and Puig, V., 2016. Stochastic model predictive control based on Gaussian processes applied to drinking water networks. IET Control Theory & Applications, 10(8), pp.947‐955. Wang, C.Y., Doan, M.L., Xue, L. and Barbour, A.J., 2018. Tidal response of groundwater in a leaky aquifer—Application to Oklahoma. Water Resources Research, 54(10), pp.8019‐8033. Walck, C., 2007. Hand‐book on Statistical Distributions for Experimentalists, University of Stockholm Internal Report SUF‐PFY/96‐01. 103 Yeh, H.D. and Huang, Y.C., 2005. Parameter estimation for leaky aquifers using the extended Kalman filter, and considering model and data measurement uncertainties. Journal of Hydrology, 302(1‐4), pp.28‐45. Yeh, T.C.J. and Liu, S., 2000. Hydraulic tomography: Development of a new aquifer test method. Water Resources Research, 36(8), pp.2095‐2105. 104