Modelling the vertical flow of groundwater in a saturated-unsaturated medium

(1)

MODELLING THE VERTICAL FLOW OF

GROUNDWATER IN A SATURATED-UNSATURATED

MEDIUM

RENDANI VELE MAKAHANE

STUDENT NO.: 2012033987

Submitted in fulfilment of the requirements for the degree

Doctor of Philosophy majoring in Geohydrology

in the

Faculty of Natural and Agricultural Sciences

(Institute for Groundwater Studies)

at the

University of the Free State

Supervisor: Prof. Dr. Abdon Atangana

BLOEMFONTEIN

January 2020

(2)

i

DECLARATION

I, Rendani Vele Makahane, hereby declare that the dissertation hereby submitted by me to the Institute for Groundwater Studies in the Faculty of Natural and Agricultural Sciences at the University of the Free State, in fulfilment of the degree Doctor of Philosophy, is my own independent work and I have not previously submitted it for a qualification at another institution of higher education. In addition, I declare that all sources cited have been acknowledged by means of a list of references.

I furthermore cede copyright of the dissertation and its contents in favour of the University of the Free State.

In addition, the following papers have been submitted to Springers:

1. Makahane R.V., & Atangana A., 2020. Analysis of the existing model for the vertical flow of groundwater in the saturated-unsaturated zone.

2. Makahane R.V., & Atangana A., 2020. New model of the saturated-unsaturated groundwater flow with power law and scale-invariant mean square displacement.

3. Makahane R.V., & Atangana A., 2020. New model of the 1-d unsaturated-saturated groundwater flow with crossover from usual to confined flow mean square displacement. 4. Makahane R.V., & Atangana A., 2020. New model of the 1-d unsaturated-saturated

groundwater flow with crossover from usual to sub- flow mean square displacement. 5. Makahane R.V., & Atangana A., 2020. New model of the 1-d saturated-unsaturated

groundwater flow using the fractal fractional derivative.

6. Makahane R.V., & Atangana A., 2020. Application of the fractional-stochastic approach to the saturated-unsaturated zone model.

7. Makahane R.V., & Atangana A., 2020. Transfer function of the Sumudu, Laplace transforms and their application to groundwater.

Rendani Vele Makahane January 2020

(3)

ii

DEDICATION

A special dedication to my nephews and nieces, may this motivate you to do it a little extra when it comes to educating yourself. To my late nephew Zwothe Makahane, this is yours baby.

(4)

iii

ACKNOWLEDGMENTS

It has been a great pleasure and a valuable experience to do my Ph.D. at the Institute of Groundwater Studies, University of the Free State; being assisted by a very helpful and experienced supervisor Prof. A. Atangana who was dedicated to helping during this research. May God continue to give our team wisdom and strength to learn more and share our experience with other scientists who are interested in this topic. This has created for me a platform for much more research.

Firstly, I would like to thank my supervisor Prof. A. Atangana for his advice, guidance, comments, insight, and encouragement over the course of this research. I will also like to thank the National Research Foundation for funding the 1st year of this degree and Council for Geoscience for funding the final year of this degree.

A huge thank you to Rudzani Makahane and his family, Andani Makahane and his family, and Mbavhalelo Makahane for their love, support, encouragement, and belief in me, but most importantly I want to thank my parents Mr. Tshitswiriri Richard Makahane and the late Mrs. Thidziambi Rosemary Makahane for instilling in me the values of hard-work, enthusiasm, and professionalism that I have carried with me in whatever I do.

(5)

iv

ABSTRACT

The classical equation developed to replicate the vertical groundwater flow through the saturated-unsaturated zone is based on the concept of rate of change. While a lot of success has been achieved using this concept, researchers discovered that only classical mechanical homogenous problems with no memory could be depicted using these mathematical tools. Therefore, a need to formulate new mathematical tools that can be used to depict non-conventional behaviour was urgent as some physical problems were not understandable. For this reason, we introduced to the saturated-unsaturated groundwater flow equation differential operators that are based on the power-law kernel, exponential decay kernel and Mittag-Leffler function to develop new models. These operators are the Caputo, Riemann-Liouville, Caputo-Fabrizio and the Atangana-Baleanu. Some of these operators are local and unable to model problems that display fractal behaviours. To account for this, we used the recently introduced fractal-fractional approach to develop more models. However, the above differential operators are based on sets of constant parameters/coefficients and largely ignore random fluctuations that are likely to occur in nature. To account for randomness, a stochastic approach is used. However, this approach is unable to capture fading memory and long-range behaviour. To help us capture problems with randomness, fading memory and long-range behavior we used the recently established stochastic-fractional approach to generate more models that give us a better understanding of the complex flow. The numerical solutions and stability analysis for all the models are detailed. We also suggested a new approach for modelling groundwater flow based on the Sumudu and Laplace transform. The new models presented in this thesis will be relevant to the majority if not all problems associated with unsaturated-saturated groundwater flow. For this reason, we have not focused on any specific problem/study area, however, theoretical parameters and assumptions are presented to illustrate the effectiveness of the new models, in which readers can make possible modifications to parameters and assumptions to match their specific groundwater problem/study area.

Keywords: Groundwater, saturated-unsaturated zone, Caputo, Caputo-Fabrizio, Atangana-Baleanu, fractal-fractional, stochastic-fractional, transfer function, Sumudu transform, Laplace transform.

(6)

v

LIST OF GREEK NOTATIONS

𝛼 Alpha 𝜔 Omega

𝛽 Beta 𝜋 Pi

Δ Delta 𝜓 Psi

𝜖 Epsilon variant 𝜏 Tau

𝛤 Gamma 𝜃 Theta

Λ Lambda 𝜎 Sigma

LIST OF ABBREVIATIONS AND NOTATIONS

𝜓_𝑏 Air entry suction pressure head

1-d One dimension

𝐴𝐵 Atangana-Baleanu 𝜃₀ residual water

content

𝐴𝐵𝐶

Atangana-Baleanu-Caputo

𝐾_𝑠 saturated hydraulic conductivity

𝐶𝐹 Caputo-Fabrizio 𝜃_𝑠 saturated water

content

𝐹𝐹𝑅 _{Fractal Fractional}

Riemann Liouville

𝜃 soil water content

𝑓 Function 𝑆_𝑠 specific storage of

the soil 𝜆 pore size distribution

index

𝑡 Time

(7)

vi

List of figures

Figure 1: Schematic representation of the research framework ... 7 Figure 2: Integral transforms (La Rosa, n.d) ... 10 Figure 3: Representation of a transfer function (HELM, 2008) ... 11 Figure 4: Relationship between the pressure head and volumetric water content for

saturated-unsaturated flow (modified after Nishigaki and Kono, 1980)... 17 Figure 5: Comparison between computational molecules for the 3 schemes (Recktenwald, 2011). ... 30 Figure 6: Transfer function of the model with the classical equation using the Laplace transform in time indicates that the system behaves like a high pass filter ... 142 Figure 7: Transfer function of the model with the classical equation using the Laplace transform in space indicates that the system behaves like a low pass filter ... 142 Figure 8: Transfer function of the model with the Caputo derivative using the Laplace transform shows that the system behaves like a low pass filter. ... 143 Figure 9: Transfer function of the model with Caputo-Fabrizio derivative using the Laplace transform shows that the system behaves like a low pass filter. ... 143 Figure 10: Transfer function of the model with the classical equation using the Sumudu transform in space indicates that the system behaves like a high pass filter ... 142 Figure 11: Transfer function of the model with the classical equation using the Sumudu transform in space indicates that the system behaves like a low pass filter ... 142 Figure 12: Transfer function of the model with Caputo derivative using the Sumudu transform with the indicates that the system behaves like a high pass filter. ... 143 Figure 13: Transfer function of the model with the Sumudu transform using the Caputo-Fabrizio

derivative shows that the system behaves like a high pass filter. ... 143 Figure 14: Transfer function of the model with Atangana-Baleanu derivative using the Laplace transform indicates that the system behaves like a low pass filter. ... 144 Figure 15: Transfer function of the model with Atangana-Baleanu derivative using the Sumudu transform indicates that the system behaves like a high pass filter. ... 144 Figure 16: Numerical simulations for the classical saturated groundwater flow equation... 147 Figure 17: Numerical solution of the saturated groundwater flow equation using fractal and power law kernel with 𝜷 = 𝟎. 𝟗 and 𝜶 = 𝟏 ... 148 Figure 18: Numerical solution of the saturated groundwater flow equation using fractal and the power law kernel with 𝜷 = 𝟏 and 𝜶 = 𝟎. 𝟗 ... 149 Figure 19: Numerical solution of the saturated groundwater flow equation using fractal and power law kernel with 𝜷 = 𝟎. 𝟗𝟓 and 𝜶 = 𝟏 ... 150 Figure 20: Numerical solution of the saturated groundwater flow equation using fractal and power law kernel with 𝜷 = 𝟎. 𝟗 and 𝜶 = 𝟎. 𝟗 ... 151 Figure 21: Numerical solution of the saturated groundwater flow equation using fractal and power law kernel with 𝜷 = 𝟎. 𝟗𝟓 and 𝜶 = 𝟎. 𝟗𝟓 ... 152 Figure 22: Numerical solution of the saturated groundwater flow equation using fractal and power law kernel with 𝜷 = 𝟏 and 𝜶 = 𝟏 ... 153

(12)

xi

Figure 23: Numerical solution of the saturated groundwater flow equation using fractal and power law kernel with 𝜷 = 𝟎. 𝟒𝟓 and 𝜶 = 𝟏 ... 154

(13)

1

CHAPTER 1 INTRODUCTION

1.1. BACKGROUND STUDY

With the rapid expansion of the population, agriculture, and industries, there is a substantial increase in the usage of groundwater resources to supplement surface water supplies. Due to this profound effect, groundwater development should be undertaken only after careful planning, which includes analyzing, understanding, and modeling the groundwater movement. Therefore, it is desired to have a model that predicts the flow of water through the saturated zone in both confined and unconfined aquifers and for circulation in the unsaturated zone. Instead of dealing with the complexity of the real world, certain assumptions can be introduced in the form of a mathematical model consisting of equations which connect parameters and variables relevant to the physical process (Pelka, 1983). For this study, the physical process under consideration takes place in the geological formation above an aquifer, and in an aquifer, this area is referred to as the saturated-unsaturated zone. An aquifer is anatural system and has a high degree of variability (Peter et al., 2011). Based on its structure, hydraulic performance, texture, geological formation and the mobility of water, three main types of aquifers are identified; unconfined, confined and semi-confined/leaky (Christiansen &Hamblin, 2014).A confined aquifer has an aquitard at its top and its base; the groundwater here is under pressure. An unconfined aquifer consists of an unsaturated medium and a saturated medium separated by a water table. In a semi-confined aquifer, one boundary is an aquiclude, and the other is an aquitard or both boundaries are aquitards, (Gushman & Tartakovsky, 2016; Pelka, 1983). Aquifers occur in a myriad of the geological formation of either unconsolidated and consolidated bedrock or both. Unconsolidated geological materials are made of porous media such as gravel, sand, silts, and clay; in this type of formation, groundwater is stored and travels in pore spaces between the particles. In consolidated bedrock formations, the water is stored and travels in fractures, fissures, and joints in the rocks, this is called secondary porosity (Uhl et al., 2009). The process of groundwater movement is a difficult one due to the presence of two distinct regions, namely, unsaturated and saturated zone. In the next few paragraphs, the essentials of the theory of subsurface water flow in the saturated and unsaturated area are discussed.

(14)

2

1.1.1 Theory of groundwater flow

The vadose zone, also referred to as the unsaturated zone is found immediately below the surface and above the groundwater level, and it is essentially a two-phase flow of two immiscible fluids containing both water and air in the pores (Freeze, 1971; Kresic, 2007). The unsaturated zone is further divided with respect to the occurrence and circulation of water into the uppermost zone of soil water, gravitational water, and the capillary fringe (Kresic, 2007). Flow in this zone is determined by the negative hydrostatic potential and the gravitational potential owing to capillary suction (Allepalli & Govindaraju, 2009). As water seeps into the ground, it substitutes the soil moisture deficit, and the excess water moves down by the force of gravity, which then builds up the groundwater (Gushman & Tartakovsky, 2016). Below the water table, groundwater moves both vertically and laterally, and water constantly occupies all pore spaces and/or rock fractures; this zone is referred to as the groundwater zone or saturated zone. In the saturated zone, movement is overseen by the total energy head or the piezometric head, made up of the positive pressure potentials and the elevation head (Allepalli & Govindaraju, 2009).

The phenomena by which water is absorbed or released from storage by the soil differ significantly in the saturated zone and in the unsaturated zone. In the unsaturated zone, water is released from storage primarily through a process called desaturation, although some distortion of the soil’s framework may occur. In the saturated area, water is taken into or released from storage because of the distortion of the soil’s framework and a consequent change in the porosity of the soil. In dealing with the unsaturated flow, soil physicists generally tend to ignore the deformation characteristics of the soil skeleton, while the soil engineers and hydrogeologists concentrating on the saturated zone, frequently neglect the effects of desaturation (Narasimhan, 1975). The water in the saturated area is enough to be withdrawn for drinking, irrigation or other uses; therefore, this zone is referred to as an aquifer (Mahajan, 2008).

When solving a natural problem, for instance, the movement of water through a porous media, knowledge on aquifer properties is essential (Fitts, 1990). Aquifer properties that govern the groundwater flow include porosity (n), permeability, storativity (S), hydraulic conductivity (K), hydraulic head (h), hydraulic gradient (i), transmissivity (T) and specific yield (Sy). Porosity is an essential property of rocks that allow the storage and flow of water underground. It directly influences the permeability, the hydraulic conductivity of rocks, therefore the velocity of

(15)

3

groundwater, and other fluids that may be present (Mahajan, 2008). When modelling the movement of groundwater, the content of soil water and the hydraulic conductivity are essential parameters. The hydraulic conductivity of a porous material is subject to the shape, connectedness, and size of the filled pores. It is shallow at low and moderate water content, and it increases non-linearly to its saturated value as water content increases (Hu et al., 2015). The hydraulic conductivity under saturated conditions is determined by the soil grain size (porosity); meanwhile, it is determined by porosity and degree of saturation in the unsaturated zone. Using the properties mentioned above, mathematical models for groundwater flow can be generated.

1.1.2 Equation of motion

The fundamental equation of groundwater flow developed initially for only saturated conditions was first established in 1856 by Henri Darcy (Hsu et al., 2002). Laplace’s equation governs the Darcian macroscopic theory for the hydraulic head (Tritscher et al., 2001). Darcy’s law states that for a given type of porous medium, the rate of flow q is proportional to the change in head ℎ and inversely proportional to the distance of 𝑧, calling the proportionality constant K the hydraulic conductivity resulting to the following equation (Wang & Anderson, 1982):

𝑞 = −𝐾∂ℎ

∂z (1.1) Darcy’s law is applied to most groundwater flow. However, it is not applicable if the medium is too irregular or if the flow velocity is too high in a medium with large pores (Fitts, 1990).

1.1.3 Equation of continuity

The conservation law, as applied to fluid transfer, is also known as an equation of continuity. This law states that under transient flow, the rate at which saturation changes in a closed system is the same as the rate of change of the total amount of fluxes in and out of the system (Pelka, 1983). The equation of continuity that describes the transient movement of water in a vertical cross-section of a saturated-unsaturated aquifer under isothermal conditions and in the absence of sources and sinks is given by (Shahraiyni and Ashtiani, 2008; Hsu et al., 2002):

𝜕𝜃

(16)

4

where ∇ is the gradient operator 𝜕 𝜕𝑧⁄ ; 𝑧 is the vertical direction. 𝜃 is the volumetric content of soil-water (L0_{) defined as 𝑛𝑆}

𝑤, 𝑛 is the medium’s porosity (L0), 𝑆𝑤 is the level of volumetric fluid

saturation, 0 ≤ 𝑆_𝑤 ≤ 1 measured in (L0_{); 𝑞 is specific discharge (LT}-1_);_{𝑡 is time (T). Richards} was perhaps the first person to express the equation of continuity for transient soil-water flow in a formula of a partial differential, parabolic equation in the soil physics literature. Familiarly known as Richard’s equation, it can be written as:

𝜕𝜃 𝜕𝑡 = 𝜕 𝜕𝑧[𝐾(𝜓) ( 𝜕𝜓 𝜕𝑧 + 1)] (1.3) where 𝐾 represents the hydraulic conductivity, 𝜃 represents the volumetric water content, 𝜓 represents the pressure head, 𝑧 represents the vertical direction, and time is represented by 𝑡. The Richards’ equation numerical solution has a disadvantage of being computationally expensive and unpredictable since it is not guaranteed that a solver will converge for a certain set of constitutive soil relations (Tocci et al., 1997). Therefore, this method cannot be applied where there is a high risk of non-convergence. Another disadvantage is over-emphasising of the role of capillary and being overly simplistic (Gray & Hassanizadeh, 1991).

Water movement through a porous media has been of great interest to humankind since early history (Nassehzadeh-Tabrizi et al., 1978). The study of water flow in the vadose /unsaturated zone and the saturated zone has frequently been accomplished by separate equations, applied distinctly on each zone. In most cases, the saturated zone is found in the area of interest, and it can fluctuate depending on how the water table moves between the unsaturated and saturated zone. In this case, Richards’ equation alone cannot govern the flow in both zones (Cheng & Gulliksson, 2003). To model a flow region where both saturated and unsaturated regimes coexist, one should combine the physical features governing the nature of each of these crucial phenomena (Narasimhan, 1975). Up until recently, very little attention seems to have been paid to the expansion of the combined saturated-unsaturated model, considering the factors particular to each regime of flow.

1.1.4 Principal equation of flow in porous media

Most of the groundwater flow models cover only either the saturated or unsaturated part of an aquifer, treating both parts separately. The flow of groundwater through the saturated-unsaturated zone is tightly coupled, and this coupling is vital in addressing problems of flow. It becomes

(17)

5

evident that groundwater flow in the saturated and unsaturated must be viewed as a whole and should be described by a sole mathematical model. By considering the assumptions that air pressure remains constant, i.e. zero in the unsaturated zone and there is negligible flow of water vapour, the flow of water through the saturated-unsaturated media can be mathematically represented by the combination the generalized Darcy’s Law and the continuity equation (List and Radu, 2015; Pelka, 1983;Zimmerman & Bodvarsson, 1989). The equation includes the hydraulic properties of the soil which are a function of the suction head of the soil and therefore is non-linear (Allepalli and Govindaraju, 1996). The 1-d saturated-unsaturated groundwater flow equation is given below: [𝑆_𝑆𝑆_𝑎(𝜓) + 𝐶(𝜓)]𝜕𝜓 𝜕𝑡 = 𝜕 𝜕𝑧[𝐾𝑧(𝜓) ( 𝜕𝜓 𝜕𝑧 − 1)] (1.4) where: 𝜓 = pressure head, 𝑛= porosity, 𝑆_𝑆= specific storage of the soil, 𝑆_𝑎 = saturation of the aqueous phase= 𝜃

𝑛 , 𝐶(𝜓) = capillary capacity of soil = 𝑑𝜃 𝑑𝜓⁄ , 𝑞(𝑡) = groundwater flux

=−𝐾𝑧(𝜓) ( 𝜕𝜓

𝜕𝑧− 1), 𝐾𝑧= hydraulic conductivity, 𝑧 = vertical co-ordinate.

1.2. PROBLEM STATEMENT

The mathematical model in equation (1.4) was constructed to replicate the vertical groundwater flow through the saturated-unsaturated zone. The geological structures of a saturated zone are different from those of unsaturated zone thus; the groundwater flows from one geological formation to another with different scales. In several already published materials, it was argued that the differential operator based on the concept of rate of change applied in classical mechanics by Sir Newton could not accurately replicate the movement of materials in different scales. This is due to the fact that the differential operator is local, thus cannot reproduce the long-range behavior of a flow in a fracture or cannot describe the flow from one scale to another. It is also confirmed that this mathematical tool “derivative (rate of change)” does not produce a mean square displacement with crossover behavior. In this light, it is strongly argued that the above mathematical model will not be able to replicate the vertical movement of groundwater flow in the saturated-unsaturated zone. We strongly support this argument due to the fact that the geological formation within which this process takes place is not homogeneous, and the above model does not in any way incorporate into its mathematical formulation the effect of heterogeneity, neither the scaling factor of the aquifer. For this study, we will use the newly established differential

(18)

6

operators based on the exponential decay kernel, power-law kernel and the Mittag-Leffler function to develop new models. New models will also be derived using the recently introduced fractal-fractional and the fractal-fractional-stochastic approach. A new approach in groundwater modelling using the Sumudu and Laplace transform is also suggested.

1.3. AIMS AND OBJECTIVES

1.3.1 Aims

This study aims to model the vertical flow of groundwater in a saturated-unsaturated medium. 1.3.2 Objectives

 Analyse existing models for vertical flow in saturated-unsaturated zone

 Develop a new model with power-law and fractal properties

 Develop a new model with crossover from usual to confined flow mean square displacement

 Develop a new model with crossover from usual to sub-flow mean square displacement

 Apply the fractional-stochastic approach to the 1-d saturated-unsaturated zone model

 Sumudu, Laplace transforms and their application to groundwater

 Numerical simulations and interpretation

1.4. RESEARCH OUTLINE

 Chapter 1: Introduction

 Chapter 2: Literature review

 Chapter 3: Analysis of an existing model of vertical flow in saturated-unsaturated zone

 Chapter 4: New model with power-law and fractal properties

 Chapter 5: New model with crossover from usual to confined flow mean square displacement

 Chapter 6: New model with crossover from usual to sub-flow mean square displacement

 Chapter 7: Application of the fractional-stochastic approach to saturated-unsaturated zone model

 Chapter 8: Sumudu, Laplace transforms and their application to groundwater

(19)

7

1.5. RESEARCH FRAMEWORK

Figure 1: Schematic representation of the research framework

Modelling the vertical flow of water in an unsaturated-saturated medium

Flow of water in the unsaturated zone

Mathematical model for water flow in unsaturated-saturated zone

Analysis of existing models of vertical flow in unsaturated-saturated zone

Darcy’s Law

Flow of water in the saturated zone

Richards’s Law

Application of the stochastic-fractional

approach to

unsaturated-saturated zonemodel

New model with power-law and fractal properties

New model with crossover from usual

to confined flow mean square displacement

New model with crossover from usual

to sub-flow mean square displacement

Sumudu, Laplace transforms and their application to groundwater

Developing new models

Numerical simulations, discussion, results & conclusion

(20)

8

CHAPTER 2 LITERATURE REVIEW ON FRACTIONAL CALCULUS AND STOCHASTIC

APPROACH

2.1. INTRODUCTION

To capture heterogeneities observed by mankind in their day-to-day life, mathematicians count on mathematical models to produce analytical and numerical solutions. To solve equations analytically, integral transforms are used; meanwhile, fractional derivatives are used to solve equations numerically. Usually, these models are built with constant coefficients; however, it was discovered by statisticians that changing a constant coefficient into a distribution assistance with capturing processes with statistical settings (Atangana & Araz, 2019). This section serves as a brief propaedeutic chapter for understanding the methods, results, and conclusions presented in this thesis.

2.2. FRACTIONAL CALCULUS

History indicates that Isaac Newton and Gottfried Wilhelm Leibniz discovered calculus autonomously in the 17th century. In his discovery of calculus, G. W. Leibniz introduced the idea of a symbolic technique and used the symbol 𝑑𝑛_{𝑦 𝑑𝑥}_⁄ 𝑛 _{= 𝐷}𝑛_{𝑦 for the n}th_{derivative, where n is a}

non-negative integer (Debnath, 2004). In 1695, Guillaume de L’Hospital wrote to G. W. Leibniz and asked him what will the result be if 𝑛 = 1/2. Leibniz responded using the following words “An apparent paradox, which one-day useful consequences will be drawn.” Those are the words that marked the beginning of fractional calculus (Loverro, 2004; Nishimoto, 1991). Fractional calculus is an arena of mathematics that raises out of the traditional descriptions of the calculus derivative and also integral operators of non-integer order (Petras, 2011; Dalir & Bashour, 2010). It can be used to model engineering and physical processes, which are described best using fractional differential equations. However, it is important to note that in many cases, standard mathematical models of integer-order derivative, together with non-linear models, do not work efficiently. Lately, fractional calculus has been used in many other disciplines such as economics, chemistry, continuum mechanics, physics, biology, electromagnetics, image processing etc. (Atangana & Secer, 2013; Dalir & Bashour, 2010). Researchers such as Laplace, Fourier, Euler are one of the few that experimented with fractional calculus and mathematical outcomes. While using their methodology and notation, many found definitions that are appropriate for the concept of a non-integer order derivative or integral (Loverro, 2004; Shukla & Sapra, 2019). Such

(21)

9

definitions include integrals and derivatives of, Caputo, Riemann-Liouville, Caputo-Fabrizio, and the recently introduced Atangana-Baleanu. The fractional calculus can be understood precisely by knowing some simple mathematical definitions like Mittag- Leffler function, Laplace transforms, Beta function and Gamma function. These are discussed in the following subsections.

2.2.1 Gamma function

The gamma function (Γ) appears to be a crucial function of fractional calculus and is the extension of the factorial function for the non-negative integers with its argument shifted down by 1, to real and complex numbers (Loverro, 2004; Shukla & Sapra, 2019). The gamma function can be represented as:

Γ(𝑧) = (𝑧 − 1)! (2.1) If the real part of the complex number 𝑧 is greater than zero, then the integral is given as:

Γ(𝑧) = ∫ 𝑥𝑧−1_𝑒−𝑥_𝑑𝑥 ∞

0

(2.2) 2.2.2 Beta function

Euler and Legendre are the first people that studied the Beta function/Euler integral. For all real number 𝑥 > 0 and 𝑦 > 0 the Beta function is defined as (Loverro, 2004):

𝛽(𝑦, 𝑧) = ∫ 𝑡𝑦−1(1 − 𝑡)𝑧−1 1

0

(2.3)

The main property of this function is its relationship to the gamma function:

𝛽(𝑦, 𝑧) = Γ(𝑦)Γ(𝑧)

Γ(𝑦 + 𝑧) (2.3.1) The above relationship is significant for fractional calculus (Shukla & Sapro, 2019).

2.2.3 Mittag-Leffler function

This function depends on complex parameters 𝛼 and 𝛽. It was named after Gosta Mittag-Leffler. If the real parts of 𝛼 are strictly positive, then the function is defined by the series below (Kowankar & Gangal, 1996; Mainardi, 2010):

(22)

10 𝐸_𝛼,𝛽(𝑧) = ∑ 𝑧 𝑘 Γ(𝛼𝑘 + 𝛽) ∞ 𝑘=0 (2.4)

In the concept of integer-order differential equations, the exponential function 𝑒𝑧_{plays an}

imperative part, such as finding the solution to 𝑦1_{+ 𝑦 = 0. In the same way, the above}

Mittag-Leffler function is very vital for the fractional-order calculus, and it has an important role in finding solutions to the non-integer order differential equations (Shukla & Sapra, 2019).

2.2.4 Laplace Transform

Pierre-Simon Laplace is an astronomer and mathematician who introduced the Laplace transform. The Laplace transform is a compelling method and it is used to solve convolution and differential equations. This integral transform allows one to turn complicated equations into simple ones (Luchko, 2019; Zachary & Tseng, 2008). Suppose you have differential equation 𝑔(𝑡), if the Laplace transform is used on an equation, it will convert it into an algebraic one 𝐺(𝑠). Using the inverse Laplace transform, the solution 𝑌(𝑠) can be mapped back to the original domain 𝑦(𝑡) (Bryant, 2008; Sontakke & Shaikh, 2015). This procedure is represented by the diagram below (La Rosa, n.d):

There are two main advantages of the Laplace transform. Firstly, problems are solved more directly; there is no need for finding a general solution of an initial value problem, and nonhomogeneous ordinary differential equations are not required to convert into homogeneous ODEs. Secondly, the use of a unit step function and Dirac’s delta makes the method prevailing for

Transform Algebraic equation 𝐺(𝑠) Differential equation g(𝑡) Inverse transform 𝑦(𝑡) 𝑌(𝑠) Difficult solution Easiest solution

(23)

11

complications with inputs (Shukla & Sapra, 2019). The formal definitions of the Laplace transform and convolution are presented in the following chapter.

2.2.4.1 Transfer function and block diagram

Another way in which the Laplace transform is used in modelling complex processes is by making use of a transfer function and representing it using a block diagram. A transfer function of a time-invariant, linear system is a proportion of the output variable (response function), 𝑍(𝑠) = ℒ{𝑧(𝑡)} to the input variable also referred to as adriving function 𝑌(𝑠) = ℒ{𝑦(𝑡)}, of the Laplace transform, with all zero initial conditions (Zwart, 2004). The transfer functions together with Bode diagrams, it can be used to describe the relationship between the cause and effect throughout a system (Aziz, 2010). By taking the Laplace transform of a differential equation, the equation is then transformed into:

𝐺(𝑠) ≡𝑍(𝑠)

𝑌(𝑠) (2.5) where, 𝐺(𝑠) is a transfer function. Linking the output variable to the input variable can be designated schematically by the use of a block diagram, as shown below:

2.2.5 Definitions, advantages, and disadvantages of the fractional derivatives and integrals

Nowadays, the fractional derivatives that are used are built based on three types of dominant calculus; the exponential decay law, Mittag-Leffler function, and power law kernel. The fractional-order derivative is not certainly exceptional; with that said, there are some accepted and standard definitions in the literature (Matlob & Jamali, 2017). In this section, we mention these definitions and also highlight their advantages and disadvantages. It is also crucial to note that despite all these disadvantages, the fractional derivatives and integrals have been used in various fields with a lot of success. The mathematical representation of these definitions is provided in chapters that follow.

System differential equation 𝐺(𝑠)

Input Output

𝑌(𝑠) 𝑍(𝑠)

(24)

12

The Riemann-Liouville fractional derivative is built using the power law kernel (Alqahtani, 2016). With the R-L fractional derivative, a random function doesn’t have to be constant at its origin, and it doesn’t have to be differentiable (Atangana, 2018b). This derivative has particular limitations when modelling practical events with fractional differential equations. R-L fractional derivative is a non-local operator and the derivative of its constant is not zero (Atangana & Bildik, 2013; Atangana, 2017). Additionally, if a random function at the origin is continuous, then the origin of its fractional derivative will have singularity. These disadvantages lessen the application of this derivative (Atangana & Secer, 2013).

In 1967 an alternative fractional derivative also founded on the power law was presented by Caputo. It has been pointed out by many researchers that fractional derivatives based on the power-law are powerful and appropriate mathematical tools that can accurately describe the flow in a media that is elastic (Alqahtani, 2016). The greatest advantage of a Caputo derivative is its ability to incorporate the traditional boundary and initial conditions in the construction of problems. Additionally, the derivative of its constant is zero, therefore the Caputo is much appreciated to describe natural problems compared to the Riemann-Liouville (Atangana & Secer, 2013; Alqahtani, 2016). However, the Caputo fractional derivative is local and demand sophisticated conditions of regularity for differentiability, this reduces its field application (Atangana & Bildik, 2013; Atangana, 2018b).

Caputo and Fabrizio recently suggested a new derivative in the field of fractional calculus. This derivative is built using the exponential decay law and it is collectively known as the Caputo-Fabrizio fractional derivative. Some benefits of this derivative over the Caputo and R-L are that it can be used where exponential law is noticed in nature, the derivative is free from singularity (Caputo & Fabrizio, 2015). In addition, this derivative can describe heterogeneities of materials and configurations with different scales (Atangana & Alqahtani, 2016). However, issues were pointed out contrary to this derivative. Firstly, the kernel is local, secondly the anti-derivative related to it is merely just an average of the function and it's integral (Atangana, 2018a). Therefore, the operator is rather just a filter than a fractional derivative. Nevertheless, many works with great success were achieved using this new derivative. This derivative is used in groundwater studies, diffusion models, mechanical engineering, thermal science and others (Alqahtani, 2016).

(25)

13

In 2016, Atangana and Baleanu presented a non-local derivative built using the generalized Mittag-Leffler function with no singular kernel and satisfies all problems pointed out contrary to the Caputo-Fabrizio fractional derivative (Alqahtani, 2016). This new approach can depict material heterogeneity and structures with diverse scales (Atangana, 2018a). Due to its Mittag-Leffler memory, the Atangana-Baleanu derivative is able to differentiate between dynamical systems taking place in diverse scales without steady-state (Atangana & Gomez-Aguilar, 2017). Mittag-Leffler function is a good filter compared to the power-law and exponential functions, making the Atangana–Baleanu fractional derivative a prevailing mathematical tool to model complex practical problems (Atangana, 2018a).

2.3. STOCHASTIC APPROACH

There are two types of modeling that are generally used to solve groundwater problems; Deterministic and stochastic. Deterministic models suppose that all the input parameters are known with certainty, and subsequently, each parameter has a deterministic value. Nevertheless, this hypothesis is not always true owing to heterogeneities associated with hydrogeological parameters. Hence, deterministic modelling can result in weighty errors in the forecast, and thus, the purpose of numerical modelling cannot be accomplished (Baalousha & Kongeter, 2006). A possible method to enumerate uncertainty owing to the spatial heterogeneity of underground systems is to view heterogeneity in a stochastic sense (Chang & Yeh, 2010). The stochastic approach has been used heavily in mathematical biology. Birth and death processes, for instance, were originally used to study the growth of populations over time (Bressloff, 2014). A stochastic process is a sequence of a random variable [𝑋(𝑡)]whose value undergo a change over time due to probabilistic laws. The properties of a stochastic process must be established from a single time series. It is assumed that many stochastic processes are Markovian. The property of a first-order Markov process is that the dependence of forthcoming values of the process on previous values rest on only the present values and not on previous observations or values. Therefore, the present value summaries the state of the process. As a result, the present value of the process is usually called the state. This also makes sense when referring to the state of the level of an aquifer. Another assumption is that the stochastic process is stationary, meaning that the probability dispersal of the process will not change with time. The statistics used when describing the distribution of a continuous-state stationary stochastic process are the sample variance, the sample mean and several autocorrelations. To estimate the mean and the variance the following equations are used:

(26)

14 𝜇̂_𝑥 = 𝑥̅ =1 Γ∑ 𝑋𝑡 𝑇 𝑡=0 (2.6) and 𝛿2 𝑥 = 1 Γ∑(𝑥𝑡− 𝑥̅) 2 𝑇 𝑡=0 (2.7)

The mean and the variance are both independent of time 𝑡. However, stochastic processes are not always stationary, so statistical methods that depend on the stationary assumptions do not apply and the problems become more complicated (Loucks et al., 2005).

(27)

15

CHAPTER 3 ANALYSIS OF THE EXISTING MODEL OF VERTICAL FLOW IN THE

SATURATED-UNSATURATED ZONE

3.1. INTRODUCTION

Instead of dealing with the complex physical reality, a mathematical model consisting of a system of equations that connects the parameters and variables, relevant to their physical meaning is built by abstracting the process and introducing certain assumptions and simplifications. For this study, a transient 1-d vertical flow equation will be analyzed based on the physics of water transfer between the unsaturated and saturated zones. This area is limited by the soil surface and by the lower impervious boundary of the aquifer.

3.2. MATHEMATICAL MODEL

The unsaturated and saturated zone can be regarded as only one continuum. What links the saturated and unsaturated zone is that the solutions of the equation governing the flow in the unsaturated zone determine the recharge and storage information used to solve the equation governing the saturated flow (Pikul et al., 1974).The equation which is used for both zones is a combination of the continuity law and Darcy’s law. Therefore, by implementing suitable values of parameters, equation (1.5) can represent the Richards’ equation which governs flow in the unsaturated zone and the groundwater flow equation in the saturated zone. The resulting equation for a saturated flow is linear and it is non-linear for the unsaturated flow (Nishigaki & Kono, 1980; Cheng & Gulliksson 2003). Wherever flow is saturated (𝜓 > 0 and𝑝 > 0), modelling is easier as 𝜃and 𝐾remains constant for a specified soil textural class (Allepalli & Gavindraraju, 1996; Kresic, 2006). For flow in the unsaturated zone (𝜓 ≤ 0 and𝑝 ≤ 0), 𝐾 depends strongly on 𝜃. As 𝜃 decreases, larger pores empty first, the result is not only that fewer pores are filled to conduct water, but the remaining pores have a smaller radius and higher resistance to flow. As a result, 𝐾 decreases intensely as 𝜃 decreases (Nimmo et al., 2009). Due to this complexity, empirical relationships are essential to model unsaturated conductivity and compute flow in the vadose zone (Cattaneo et al., 2015). The empirical solution scheme used for this study is that of Brook and Corey (1964).

(28)

16

3.3. BROOK AND COREY MODEL

Soil hydraulic functions refer to the hydraulic conductivity function, 𝐾(𝜓) and the soil water content function, 𝜃(𝜓) that are required to explain the movement of water. Numerous functions have been suggested to empirically define the soil hydraulic properties. One of the popular models is the equations of Brook and Corey (Allepalli & Gavindraraju, 1996):

𝐾(𝜓) = 𝐾_𝑠(𝜓𝑏 𝜓) 2+3𝜆 𝜓 ≤ 0 (3.1) 𝐾(𝜓) = 𝐾_𝑠 𝜓 > 0 (3.2) 𝜃(𝜓) = 𝜃₀+ (𝜃_𝑠− 𝜃₀) (𝜓𝑏 𝜓) 𝜆 𝜓 ≤ 0 (3.3) where: 𝐾_𝑠= saturated hydraulic conductivity, 𝜓 = pressure head, 𝜓_𝑏= air entry suction pressure head, 𝜃(𝜓) = soil water content, 𝜃₀= residual water content, 𝜃_𝑠= saturated water content, 𝜆 = pore size distribution index.

3.4. PARAMETERS FOR SATURATED-UNSATURATED ZONE

The functional relationship between 𝜓 and 𝜃 is called the water retention curve and it is shown in figure 4 below. Incorporated in this figure are the soil hydraulic functions ( 𝐾(𝜓) and 𝜃(𝜓)) showing how they differ from unsaturated to the saturated zone. This figure is modified after Nishigaki and Kono, 1980.

(29)

17

3.5. ANALYTICAL TOOLS FOR SATURATED-UNSATURATED POROUS MEDIA

Now that the flow in the unsaturated and saturated zone has been distinguished, classical equations suitable for each zone will be derived below using powerful methods that result in explicit forms of a solution. In the case of the saturated zone, analytical and numerical solutions will be derived using integral transform, the method of separation of variables and numerical schemes respectively. The conditions under which the numerical method used converges will be derived and presented in detail. Due to the complexity of the unsaturated model, we relied only on the numerical method to derive an approximate solution. Some numerical simulations will be depicted for different values of the used input parameters. To familiarize readers with these methods, below we will discuss the definition and some properties of the Fourier Transform and Laplace Transform that will be useful for understanding this section.

3.5.1 Definition and properties of the Laplace Transform

Definition: If the Laplace transform of𝑓(𝑡), represented by 𝐹(𝑠) orℒ{𝑓(𝑡)} is defined for [0, ∞), then the Laplace integral is given by:

(Saturated zone)

Volumetric moisture content (𝜃)

P re ss ur e he ad (𝜓 ) 𝐾(𝜓) = 𝐾_𝑠 (Unsaturated zone) 𝜃 = 𝑛 ∴ 𝜃 𝑛 = 1 𝐾(𝜓) = 𝐾_𝑠(𝜓𝑏 𝜓) 2+3𝜆 𝜃(𝜓) = 𝜃₀+ (𝜃_𝑠 − 𝜃₀) (𝜓𝑏 𝜓) 𝜆 0 𝜃 𝑛= 0 𝜓 > 0, 𝑝 > 0 𝜓 ≤ 0, 𝑝 ≤ 0

Figure 4: Relationship between the pressure head and volumetric water content for saturated-unsaturated flow (modified after Nishigaki and Kono, 1980).

(30)

18

ℒ{𝑓(𝑡)} = 𝐹(𝑠) = ∫ 𝑒−𝑠𝑡_{𝑓(𝑡) 𝑑𝑡} ∞

0

(3.4)

The Laplace transform has a number of properties. We state few definitions of these properties in two ways, first in words and then in symbols just like we did for the above defined Laplace transform.

The first very useful property is the linearity of the Laplace transform:

ℒ is a linear operator. This means that of any two functions 𝑓 and 𝑔 for which the Laplace is defined and two constants 𝑎 and 𝑏 ∈ R we have:

ℒ{𝑎𝑓 + 𝑏𝑔} = 𝑎ℒ{𝑓} + 𝑏ℒ{𝑔} (3.4.1) If ℒ{𝑓} = 𝐹(𝑠) and ℒ{𝑔} = 𝐺(𝑠) then we have:

ℒ{𝑎𝑓 + 𝑏𝑔} = 𝑎ℒ{𝐹(𝑠)} + 𝑏ℒ{𝐺(𝑠)} (3.4.1.1)

Laplace transform of a first derivative. If ℒ{𝑓} = 𝐹(𝑠) then

ℒ[𝑓′(𝑡)] = 𝑠ℒ[𝑓(𝑡)] − 𝑓(0) (3.4.2)

Laplace transform of a second derivative. If ℒ{𝑓} = 𝐹(𝑠) then:

ℒ[𝑓′′(𝑡)] = 𝑠2ℒ[𝑓(𝑡)] − 𝑠𝑓(0) − 𝑓′(0) (3.4.3)

Time scaling. Let ℒ{𝑓(𝑡)} = 𝐹(𝑠) then

ℒ{𝑓(𝑎𝑡)} =1 𝑎𝐹 (

𝑠

𝑎), (3.4.4) For any𝑎 > 0

To prove this property, we start with replacing some variables in the definition of the Laplace transform; we use the change of variable 𝑎𝑡 = 𝜏, from where𝑑𝑡 = 𝑑𝜏/𝑎. Note that for 𝑎 > 0 the limits of integration will not change:

ℒ{𝑓(𝜏)} =1 𝑎∫ 𝑓(𝜏)𝑒 −_𝑎𝑠𝜏_{𝑑𝜏 =} 1 𝑎𝐹 ( 𝑠 𝑎) ∞ 0 (3.4.4.1)

(31)

19

Note: if we allow any sign for 𝑎 ≠ 0 then

ℒ{𝑓(𝑎𝑡)} = 1 |𝑎|𝐹 (

𝑠

𝑎) (3.4.4.2)

Differentiation of the frequency (multiplication by t): Let ℒ{𝑓(𝑡)} = 𝐹(𝑠). Then

ℒ{𝑡𝑓(𝑡)} = −𝐹′(𝑠)

To prove this property, we differentiate either side of the Laplace transform with respect to 𝑠 to get:

𝐹′_{(𝑠) = ∫ (−𝑡)𝑒}∞ −𝑠𝑡_{𝑓(𝑡) 𝑑𝑡}

0 (3.4.5)

Time delay: Letℒ{𝑓(𝑡)} = 𝐹(𝑠). Then

ℒ{𝑓(𝑡)} = 𝑒−𝑒𝑇_{𝐹(𝑠) (3.4.6)}

The property time delay implies that the function 𝑓(𝑡) in the Laplace transform is delayed by 𝑇 seconds and “zero padded” up to𝑇. Therefore, 𝑓(𝑡) = 𝑓(𝑡 − 𝑇):

ℒ{𝑓(𝑡)} = ∫ 𝑒−𝑠𝑡𝑓(𝑡 − 𝑇) 𝑑𝑡 ∞ 0 (3.4.6.1) Let: 𝑡 = 𝜏 + 𝑇 ℒ{𝑓(𝑡)} = ∫ 𝑒−𝑠(𝜏+𝑇)_{𝑓(𝜏) 𝑑𝑡} ∞ 0 (3.4.6.2) ℒ{𝑓(𝑡)} = ∫ 𝑒−𝑠𝑇𝐹(𝑠) ∞ 0 (3.4.6.3)

Convolution: The convolution assists in solving integral equations. If we let 𝐹(𝑠) and 𝐺(𝑠) denote

the transforms of 𝑓(𝑡) and𝑔(𝑡), respectively, then the inverse of the product 𝐹(𝑠)and 𝐺(𝑠) is given by the functionℎ(𝑡) = (𝑓 ∗ 𝑔)(𝑡). It is called the convolution of 𝑓(𝑡)and 𝑔(𝑡) and can be regarded as a generalised product of 𝑓(𝑡) and𝑔(𝑡) (Mathews & Howell, 2012):

ℎ(𝑡) = (𝑓 ∗ 𝑔)(𝑡) = ∫ 𝑓(𝜏) 𝑔(𝑡 − 𝜏)𝑑𝜏

𝑡

0

(32)

20

The integral representation of the Laplace transform of a convolution is therefore given by: Let: ℒ{𝑓(𝑡)} = 𝐹(𝑠) and ℒ{𝑔(𝑡)} = 𝐺(𝑠) then,

ℒ[𝑓(𝑡) ∗ 𝑔(𝑡); 𝑠] = ℒ[𝐹(𝑠)]ℒ[𝐺(𝑠)] (3.4.7.1) 3.5.2 Definitions and properties of the Fourier Transform

Properties of the Fourier transform include linearity, time delay, time scaling, duality, correlation, time derivative, convolution, etc. However, we will only give the definition of the Fourier Transform.

Definitions: The Fourier transform of a function 𝑓 is usually represented by adding a circumflex 𝑓̂ . There are many forms of stating the Fourier transform of an integrable function𝐹: 𝑅 → 𝐶. In our case, the following definition will be used:

𝑓̂(𝜔) = ∫ 𝑓(𝑡)

∞

−∞

𝑒−𝑖𝜔𝑡𝑑𝑡 (3.5) 3.5.3 Definition of the inverse Fourier Transform

This transform has an important role in the application of the Fourier transform. Ideally, in many situations, the Fourier transform is applied first, performing some simplification, then apply the inverse Fourier transform. The definition of the inverse Fourier is given by:

ℱ−1_{𝑓̂(𝜔) =} 1 2𝜋∫ 𝑓̂(𝜔)𝑒 𝑖𝜔𝑡 ∞ −∞ 𝑑𝜔 (3.6)

3.6. SATURATED FLOW EQUATION

For flow in the saturated zone, we now know that the volumetric water content is equal to the porosity, we can, therefore, cancel them out of the equation. The hydraulic conductivity in the saturated zone is constant and the capillary capacity becomes zero, therefore the equation for vertical groundwater flow in the saturated zone is given as (Maslouhi et al., 2009):

(𝑆_𝑆)𝜕𝜓 𝜕𝑡 = 𝜕 𝜕𝑧[𝐾𝑠( 𝜕𝜓 𝜕𝑧 − 1)] (3.7) Since the hydraulic conductivity of the saturated zone is constant, the above equation can be further simplified into:

(33)

21 𝜕𝜓 𝜕𝑡 = 𝐾_𝑠 𝑆_𝑠 𝜕2_𝜓 𝜕𝑧2 (3.8)

In the next few pages, we attempt to find the analytical solution, numerical solution, and the stability of the solutions generated for equation (3.8).

3.6.1 Analytical solution using the Integral Transform

The above equation for the saturated zone (3.8) can be solved analytically using an integral transform, for instance, the Fourier transform or Laplace, for this case we will make use of both the Laplace and Fourier transform to provide an exact solution. Thus, applying the Laplace transform on either side of equation (3.8) we obtain:

ℒ (𝜕𝜓 𝜕𝑡) = ℒ ( 𝐾_𝑠 𝑆_𝑠 𝜕2_𝜓 𝜕𝑧2) (3.9)

Applying the Laplace transform of a derivative into equation (3.9) we obtain the following with 𝑠 being considered as Laplace variable:

𝑠𝜓̃(𝑧, 𝑠) − 𝜓(𝑧, 0) = (𝐾𝑠 𝑆_𝑠

𝜕2_𝜓

𝜕𝑧2) (𝑧, 𝑠) (3.9.1)

To eliminate the space component in order to obtain an algebraic, we use the Fourier transform to the above equation:

ℱ(𝑠𝜓̃ (𝑧, 𝑠) − 𝜓(𝑧, 0) ) =𝐾𝑠 𝑆_𝑠 ℱ (

𝜕2𝜓̃ (𝑧, 𝑠)

𝜕𝑧2 ) (3.9.2)

Using the properties of the Fourier transform presented earlier, the right-hand side of the equation above can further be simplified into a linear form as follows:

𝑠𝜓̃_ℱ(𝜔, 𝑠) − 𝜓_ℱ(𝜔, 0) =𝐾𝑠 𝑆_𝑠(𝑖𝜔)

2_𝜓̃

ℱ(𝜔, 𝑠) (3.9.3)

The above equation can be factorised to obtain the equation below:

(𝑠 −𝐾𝑠 𝑆_𝑠 (𝑖𝜔)

2_{) 𝜓̃}

(34)

22

We further simplify the above equation into:

𝜓̃ℱ(𝜔, 𝑠) = 𝜓_ℱ(𝜔, 0) (𝑠 −𝐾_𝑆𝑠 𝑠 (𝑖𝜔) 2₎ (3.9.5) And then: 𝜓̃_ℱ(𝜔, 𝑠) = 𝜓ℱ(𝜔, 0) (𝑆𝑠 𝐾_𝑠𝑠 + (𝜔)2) 𝐾_𝑠 𝑆_𝑠 (3.9.6)

The above formula is the multiplication of two Fourier transform of two functions; thus, one can use the convolution theorem to obtain the inverse transform as follow:

ℱ−1_(𝜓 ℱ(𝜔, 0)) = 𝜓(𝑧, 0) (3.9.7) 𝜓(𝑧, 𝑠) = ℱ−1 ( 1 𝐾_𝑠 𝑆_𝑠(𝜔2+ (√ 𝑆_𝑠 𝐾_𝑠𝑠) 2 ) ) (3.9.8) 𝜓(𝑧, 𝑠) = ℱ−1 ( 2√_𝐾𝑆𝑠 𝑠𝑠 𝐾𝑠 𝑆_𝑠 2√ 𝑆𝑠 𝐾_𝑠𝑠 (𝜔2+ (√ 𝑆𝑠 𝐾_𝑠𝑠) 2 ) ) (3.9.9) 𝜓(𝑧, 𝑠) = 1 𝐾_𝑠 𝑆_𝑠2√ 𝑆_𝑠 𝐾_𝑠𝑠 𝑒𝑥𝑝 (−√ 𝑆_𝑠 𝐾_𝑠𝑠(𝑧)) (3.9.10)

By integrating the above equation, we obtain:

𝜓(𝑧, 𝑠) = 1 𝐾𝑠 𝑆_𝑠 2√ 𝑆𝑠 𝐾_𝑠𝑠 ∫ 𝜓(𝜆, 0) 𝑧 0 𝑒𝑥𝑝 [−√𝑆𝑠 𝐾𝑠 𝑠(𝑧 − 𝜆)] 𝑑𝜆 (3.9.11)

To obtain the exact solution of the saturated zone, one can use the inverse Laplace transform on both sides of equation (3.9.11), such that:

(35)

23

𝜓(𝑧, 𝑠) = ℒ−1(𝜓(𝑧, 𝑠)) (3.9.12) 3.6.2 Analytical solution using the method of separation of variables

Using the method of separation of variable, we can now obtain an analytical solution of the above equation (3.8). This method is used to solve a wide variety of linear and homogeneous differential equations, such as the saturated groundwater flow equation. The boundary conditions are linear and homogenous, such that:

𝜓│_𝑧=0 = 𝜓│_𝑧=𝐿 = 0 (3.10) The dependence of 𝜓 on 𝑧 and 𝑡 can be represented as a product of function 𝑧 and function 𝑡, such that:

𝜓(𝑧, 𝑡) = 𝐹(𝑧)𝐺(𝑡) (3.10.1) Substitution 𝜓 back into equation (3.8) using the product rule we obtain:

𝜕(𝐹(𝑧)𝐺(𝑡) ) 𝜕𝑡 = 𝐾_𝑠 𝑆_𝑠 𝜕2_{(𝐹(𝑧)𝐺(𝑡) )} 𝜕𝑧2 (3.10.2)

Since the left-hand side is dependent only on 𝑡 and the right-hand side is dependent only on 𝑧, we can separate the variables as follows:

𝐹(𝑧)𝑆𝑠 𝐾𝑠 𝑑𝐺(𝑡) 𝑑𝑡 = 𝐺(𝑡) 𝑑2_𝐹(𝑧) 𝑑𝑧2 (3.10.3)

Variables can be further separated to obtain the following:

𝑆_𝑠 𝐾𝑠 𝑑𝐺(𝑡) 𝑑𝑡 𝐺(𝑡) = 𝑑2_𝐹(𝑧) 𝑑𝑧2 𝐹(𝑧) (3.10.4) The RHS and the LHS of equation (3.10.4) can be a function of 𝑧 and 𝑡 respectively, if they both equate to a constant value𝛼, such that:

𝑆𝑠 𝐾_𝑠

𝑑𝐺(𝑡)

𝑑𝑡 = 𝛼𝐺(𝑡) (3.10.4.1) and:

(36)

24

𝑑

2_𝐹(𝑧)

𝑑𝑧2 = 𝛼𝐹(𝑧) (3.10.4.2)

Now we try to obtain a solution that is not zero satisfying the boundary conditions. For our first attempt, we assume that for the constant 𝛼 < 0 real numbers 𝐵, 𝐶 exist such that:

𝐹(𝑧) = 𝐵𝑒√−𝛼𝑧_{+ 𝐶𝑒}−√−𝛼𝑧_(3.10.5)

From the boundary conditions given in equation (3.10), we get:

𝐹(0) = 0 = 𝐹(𝐿) (3.10.5.1) 𝐵 = 0 = 𝐶 implies that 𝜓 is identically 0, therefore, we reject this case. In the second assumption, we assume that when 𝜆 = 0 real numbers 𝐵, 𝐶 exist such that:

𝐹(𝑧) = 𝐵(𝑧) + 𝐶 (3.10.6) The same conclusion that 𝜓 is identically 0 can be made, and again we reject this case. The last case is when𝛼 > 0. Then there exists the real number 𝐴, 𝐵, 𝐶 such that:

𝐺(𝑡) = 𝐴𝑒−𝛼

𝑆𝑠

𝐾𝑠𝑡 (3.10.7)

And

𝐹(𝑧) = 𝐵 𝑠𝑖𝑛(√𝛼𝑧) + 𝐶 𝑐𝑜𝑠(√𝛼𝑧) (3.10.8) From the above information 𝐶 = 0 and that for some positive integer𝑛:

√𝛼 = 𝑛𝜋

𝐿 (3.10.9) Therefore, a general solution can be given as:

𝜓(𝑧, 𝑡) = ∑ 𝐷_𝑛 ∞ 𝑛=1 𝑠𝑖𝑛𝑛𝜋𝑧 𝐿 𝑒𝑥𝑝 (− 𝑛2_𝜋2𝑆𝑠 𝐾_𝑠𝑡 𝐿2 ) (3.10.10)

where 𝐷𝑛 is a coefficient determined by the initial condition, given the following initial

(37)

25 𝜓│𝑡=0 = 𝑣(𝑧) (3.10.10.1) We obtain: 𝑣(𝑧) = ∑ 𝐷_𝑛 ∞ 𝑛=1 𝑠𝑖𝑛𝑛𝜋𝑧 𝐿 (3.10.10.2) Multiplying both sides with 𝑠𝑖𝑛𝑛𝜋𝑧

𝐿 and integrating over [0, 𝐿] results into:

𝐷_𝑛 =2 𝐿∫ 𝑣(𝑧) 𝐿 0 𝑠𝑖𝑛𝑛𝜋𝑧 𝐿 𝑑𝑥 (3.10.10.3) Hence, the complete solution for equation (3.8) is given by:

𝜓(𝑧, 𝑡) = ∑ ( 2 𝐿∫ 𝑣(𝑧) 𝐿 0 𝑠𝑖𝑛𝑛𝜋𝑧 𝐿 𝑑𝑥) ∞ 𝑛=1 𝑠𝑖𝑛𝑛𝜋𝑧 𝐿 𝑒𝑥𝑝 (− 𝑛2𝜋2_𝐾𝑆𝑠 𝑠𝑡 𝐿2 ) (3.10.10.4) 3.6.3 Numerical solution

The analytical solutions of our model are limited to less complexities, such as the assumption of homogeneity, isotropy, simple initial condition, and simple geological formation. However, natural systems can have a more complex geological formation, a complex initial condition and they can be heterogeneous and anisotropic. Such complexities require numerical solutions. Depending on how 𝜓𝑡 is approximated, we have three basic finite difference schemes:

Crank-Nicolson, explicit, implicit and schemes. Using these schemes, we attempt to find numerical solutions for the 1-dimension saturated groundwater flow equation (3.8).

3.6.3.1 Numerical solution using forward Euler method (FTCS)

Appling the explicit forward Euler method on the equation (3.8), we obtain the following numerical solution: 𝜓𝑖 𝑛+1_{− 𝜓} 𝑖𝑛 Δ𝑡 = 𝐾_𝑠 𝑆_𝑠 [ 𝜓_𝑖+1𝑛 _{− 2𝜓} 𝑖𝑛+ 𝜓𝑖−1𝑛 (Δ𝑧)2 ] (3.11) 3.6.3.2 Numerical solution using backward Euler method (BTCS)

We obtain another numerical solution by applying the implicit backward Euler method on equation (3.8):

(38)

26 𝜓𝑖 𝑛+1_{− 𝜓} 𝑖𝑛 Δ𝑡 = 𝐾𝑠 𝑆_𝑠 [ 𝜓_𝑖+1𝑛+1− 2𝜓_𝑖𝑛+1+ 𝜓_𝑖−1𝑛+1 (Δ𝑧)2 ] (3.12) 3.7.3.3 Numerical solution using the Crank-Nicolson method

Applying the finite difference spatial discretization on equation (3.8), the implicit Crank-Nicolson discretization is as follows:

𝜓𝑖 𝑛+1_{− 𝜓} 𝑖 𝑛 Δ𝑡 = 𝐾_𝑠 𝑆_𝑠 [ 𝜓_𝑖+1𝑛+1_{− 2𝜓} 𝑖𝑛+1+ 𝜓𝑖−1𝑛+1 2(Δ𝑧)2 + 𝜓_𝑖+1𝑛 _{− 2𝜓} 𝑖 𝑛 _{+ 𝜓} 𝑖−1𝑛 (2Δ𝑧)2 ] (3.13)

3.6.4 Numerical stability analysis

In physical problems such as groundwater flow, stability analysis of an equation is very important. A procedure called von Neumann stability analysis based on the Fourier series is used to analyse the stability of finite difference schemes. The stability of finite difference schemes is linked to numerical errors. A scheme is said to be von Neumann stable if its amplification factor is less or equal to 1. Accuracy, however, requires that the amplification factor be as close to 1 as possible. In this section, the numerical solutions provided in equations 3.11-3.13 for the 1-d saturated groundwater flow equation (3.8) are subjected under von Neumann stability analysis. A detailed procedure for each solution is given below.

3.7.4.1. Stability analysis of a forward Euler method (FTCS)

Let us recall on the forward Euler method (3.11) which can also be written as:

𝜓_𝑖𝑛+1 = 𝜓_𝑖𝑛 (1 − 2𝛼) + 𝛼(𝜓_𝑖+1𝑛 + 𝜓_𝑖−1𝑛 ) (3.14) where 𝛼 =𝐾𝑠

𝑆𝑠 (

Δ𝑡 (Δ𝑧)2)

We consider a harmonic initial perturbation:

𝜓_𝑖0 _{= 𝑒}𝑖𝑘𝑚𝑧_(3.14.1)

Which with time evolves as:

𝜓_𝑖𝑛 _{= 𝜎}𝑛_𝑒𝑖𝑘𝑚𝑧_(3.14.1.1)

To analysis the stability of this scheme, let us find the amplification factor 𝜎by inserting equation (3.14.1.1) into equation (3.14):

(39)

27

𝜎𝑛+1𝑒𝑖𝑘𝑚Δ𝑧 _{= 𝜎}𝑛_𝑒𝑖𝑘𝑚𝑧(1 − 2𝛼) + 𝛼(𝜎𝑛_𝑒𝑖𝑘𝑚(z+Δ𝑧)_{+ 𝜎}𝑛_𝑒𝑖𝑘𝑚(z−Δ𝑧)) (3.14.2)

We can simplify equation (3.14.2) by pulling out the common factor, such that:

𝜎𝑛+1_𝑒𝑖𝑘𝑚Δ𝑧 _{= 𝜎}𝑛_𝑒𝑖𝑘𝑚𝑧_{[1 − 2𝛼 + 𝛼(𝑒}𝑖𝑘𝑚Δ𝑧_{+ 𝑒}−𝑖𝑘𝑚Δ𝑧_{)] (3.14.3)}

Further simplification of the above equation results into the following:

𝜎𝑛+1𝑒𝑖𝑘𝑚Δ𝑧 _{= 𝜎}𝑛_𝑒𝑖𝑘𝑚𝑧_{[1 + 2𝛼 (}𝑒𝑖𝑘𝑚Δ𝑧+𝑒−𝑖𝑘𝑚Δ𝑧

2 − 1)] (3.14.4)

Considering the definition for a hyperbolic cosine, we can rewrite the above equation (3.14.4) as: 𝜎𝑛+1𝑒𝑖𝑘𝑚Δ𝑧 _{= 𝜎}𝑛_𝑒𝑖𝑘𝑚𝑧[1 + 2𝛼(cos (𝑘

𝑚Δ𝑧) − 1)] (3.14.5)

Using the double angle identity for𝑐𝑜𝑠, equation (3.14.5) is written such that:

𝜎𝑛+1_𝑒𝑖𝑘𝑚Δ𝑧 _{= 𝜎}𝑛_𝑒𝑖𝑘𝑚Δ𝑧_{[1 + 2𝛼 (1 − 2𝑠𝑖𝑛}2₍𝑘𝑚Δ𝑧

2 ) − 1)] (3.14.6) 𝜎𝑛+1_𝑒𝑖𝑘𝑚Δ𝑧 _{= 𝜎}𝑛_𝑒𝑖𝑘𝑚Δ𝑧_{[1 − 4𝛼𝑠𝑖𝑛}2𝑘𝑚Δ𝑧

2 ] (3.14.6.1) Dividing both sides by𝜎𝑛_𝑒𝑖𝑘𝑚Δ𝑧_{we get:}

𝜎 = 1 − 4𝛼 𝑠𝑖𝑛2𝑘𝑚∆𝑧

2 (3.14.6.2) The above equation gives us the amplification factor𝜎 for the forward Euler method, and the stability condition is written as:

|1 − 4𝛼 𝑠𝑖𝑛2𝑘𝑚∆𝑧

2 | ≤ 1 (3.14.7) However, if 𝑘_𝑚∆𝑧 = 𝜋 ⇒ 𝜎 = 1 − 4𝛼. The stability of the above schema (3.14) for all 𝑘_𝑚∆𝑧 will only hold if𝛼 ≤1

2. Therefore, the forward Euler method (3.14) is conditionally stable or von

Modelling the vertical flow of groundwater in a saturated-unsaturated medium