EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS Vol. 16, No. 4, 2023, 2169-2186 ISSN 1307-5543 – ejpam.com Published by New York Business Global Analytical solution of vertically integrated equations of a confined aquifer in two-dimension using finite volume method Cordy Jourvel Itoua-Tsele1,2,∗, Antoine Filankembo Ouassissou2 1 Département des Licences, École Nationale Supérieure Polytechnique, Université Marien Ngouabi, Brazzaville, Congo 2 Département des Sciences Exacts, École Normale Supérieure, Université Marien Ngouabi, Brazzaville, Congo Abstract. Groundwater plays an important role in feeding many families and societies whose access to water is a huge problem for several reasons, then the mathematical model is an essential part of a hydrologist’s daily life to understand the physical phenomena. In this paper, we used Dupuit’s assumption to reduce the groundwater flow equation and integrated vertically it to obtain the diffusivity equation in terms of the hydraulic head, further, the introduction of fundamental ideas in groundwater research, including vocabulary and mathematics will give readers the basis knowledge. Numerous numerical techniques can be recognized, the unstructured finite volume method was used to solve a transient diffusivity equation through a coastal confined aquifer. The constructed finite volume schemes are used to study the distribution of hydraulic heads in the model domain, to determine water table fluctuations, and to predict the drawdown phenomena. 2020 Mathematics Subject Classifications: 92-10, 35A20, 35A08, 65M08 65M08 Key Words and Phrases: Confined aquifer,Vertically integrated, Dupuit’s assumption, Ground- water, finite volume method 1. Introduction Groundwater is all water found below the ground surface, in the saturation zone, and in direct contact with the ground or subsoil, and is used in many coastal regions to meet domestic, commercial, and agricultural daily needs [3] . The pore space is referred to as unsaturated when it is close to the soil surface and typically contains a mixture of air and water. The pores are constantly saturated with water as you go further into the ground. The definition of the term water table is the level in the ground above which the pore space is unsaturated and below which it is saturated, and where water that infiltrates ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v16i4.4922 Email addresses: cordy2018@yahoo.com (C. J. Itoua-Tsele), filan 64@hotmail.fr (A. F. Ouassissou) https://www.ejpam.com 2169 © 2023 EJPAM All rights reserved. C. J. Itoua-Tsele, A. F. Ouassissou / Eur. J. Pure Appl. Math, 16 (4) (2023), 2169-2186 2170 can accumulate and form an aquifer[11]. An aquifer is a geological formation which (i) contains water and (ii) permits significant amounts of water to move through it under ordinary field conditions [2]. There are two types of aquifers: (i) Confined aquifer: is an aquifer that is bounded above and below by impermeable geologic formations. Confined aquifers generally occur at significant depth below the ground surface. (ii) Unconfined aquifer: also called a water-table aquifer or phreatic aquifer , is an aquifer which has the water table as its upper boundary. Unconfined aquifers occur near the ground surface and it is directly recharged from the ground surface above it, except where impervious layers, sometimes of limited areal extent, exist between the phreatic surface (water table) and the ground surface. In addition to these, groundwater has numerous other advantages, including greater geothermal energy generation, higher water security through improved storage, and re- silience to the effects of climate change [4]. Groundwater pollution due to overexploitation for irrigation leads to a deterioration of groundwater quality. Groundwater levels are low- ered as a result of drilling operations to produce water close to or in coastal areas, which leads to saltwater intrusion [7]. [15] developed a deep learning technique for Groundwater flow equations using wells and GW-PINN without labeled data. To test the GW-PINN’s training performance under several sampling procedures and two limitations, five scenarios were created. The projected outcomes of GW-PINN were contrasted with MODFLOW and the analytical solution. The findings show that GW-PINN has a significant capacity for capturing the hydraulic head change for both confined and unconfined aquifers. [10] proposed a mathematical model to represent the transition of groundwater flow from con- fined to unconfined aquifers, specifically the classical differential operator that is based on the rate of change is replaced by a non-conventional one including the differential operator that can represent processes following the power law to capture the memory effect, after which the numerical analysis was carried out using the Caputo fractional operator, and numerical solutions were obtained. [13] presented a model of groundwater dynamics under stationary flow and solved the developed mathematical model using the Finite Difference Method, and an application on the specific study area of the Ayamonte-Huelva aquifer in Spain was performed. [5] developed a novel numerical model in order to calculate the behavior of unsteady, one-dimensional groundwater flow using the finite volume method. Various scenarios for drainage and recession from an unconfined aquifer, as well as water table fluctuations above an inclined leaky layer due to ditch recharge with a constant and variable upper boundary condition, were performed. The computed results of the explicit and implicit procedures were in good agreement with the findings of analytical approaches and laboratory tests. [14] also created a thorough analytical approach to figuring out how different rainfall recharge rates might affect groundwater flow in an unconfined sloping aquifer. The beginning water level was parallel to the impervious bottom of a gentle slope, and the unconfined aquifer’s domain was considered to be semi-infinite with an impervious bottom base. The results of the given analytical solution were compared to C. J. Itoua-Tsele, A. F. Ouassissou / Eur. J. Pure Appl. Math, 16 (4) (2023), 2169-2186 2171 those from earlier studies, and the analytical solution’s applicability was confirmed using data from a groundwater station in the Dali District of Taichung City, Taiwan, collected in 2012 and 2013. The aforementioned study models groundwater utilizing the steady conditions in some specific study regions, while the other study used unsteady groundwa- ter in one dimension. In this paper, the hydraulic behavior in a heterogeneous isotropic confined aquifer was described by a vertically integrated equation, and the analytical so- lution was carried out using the two-dimensional unstructured finite volume approach, the developed mathematical equation will be used to describe the flow of groundwater and the significance of the results of some physical phenomena considered in this study are the drawdown and water table phenomena. Steps to building a groundwater model include: • Define the domain of study; • construct a conceptual model of the groundwater (the input field data); • method of solutions; • calibration; • use the calibrated model to make predictions; • visualize and interpret the model results. 2. Dupuit’s assumption for confined aquifer flow in isotropic porous medium In order to obtain the saturated flow equation for a confined aquifer in two-dimensional space, we introduce Dupuit’s assumption. The assumption is that the flow is essentially considered to be on a horizontal plane, say the (x, y)-plane, thus: ∂ ∂z ( Kzz ∂h ∂z ) = 0, (1) to proceed we need an expression for the vector of discharge per unit width, through the entire thickness of the aquifer, Q ′ (x, y). Here, width is a line segment in the plane, and the discharge is normal to that line. Moreover, to obtain the saturated flow equation we use the following assumptions: (A-1) the aquifer is a horizontal, inhomogeneous isotropic confined aquifer; (A-2) fluid is slightly compressible and the porous medium is isotropic; (A-3) the bottom plan of the aquifer is essentially horizontal; (A-4) the confined aquifer has elastic deformation resulting from changes in hydraulic head in the vertical direction z only; C. J. Itoua-Tsele, A. F. Ouassissou / Eur. J. Pure Appl. Math, 16 (4) (2023), 2169-2186 2172 (A-5) the hydraulic conductivity and the horizontal components of flow are uniform along the entire saturated thickness of the aquifer ; (A-6) Darcy’s law is for isotropic porous medium; (A-7) the slope of the water table is small; (A-8) sinks and sources are present. 3. Mathematical description of groundwater flow 3.1. Hydraulic head The water in groundwater systems is always in motion, namely water flows from high elevation to low elevation, which is one of the primary topics that we are interested in in this study. In order to understand flow we need to understand how to express the energy in water, the hydraulic head is a crucial term to understand in groundwater. Hydraulic head is a measurement of water pressure or energy of a body of water above a given datum, or in other words, the potential energy held inside a body of water (potential energy derived from the gravitational field’s elevation of water and pressure head derived from the fluid pressure distribution) to move (see Figure 1). It’s measured in length units. The hydraulic head h is defined as follows: h = p ρ0g + z (2) where z is the elevation of the point at which the piezometric head is being considered, above some datum level (m), p and ρ are the fluid’s pressure at measurement point (Pa, force (weight) per unit area) and mass density (kg/m3), respectively, and g is the gravity acceleration (m/s2). 3.2. Darcy’s law The flux density or specific discharge q is a fictitious macroscopic rate of a flux of water passing through a unit area. Discovered experimentally for a homogeneous isotropic porous medium by Darcy in 1856, Darcy’s law is an equation that describes the flow of a fluid through a porous medium. It is used to define the relation between specific flow q and hydraulic head h and describes the flow of a fluid through a porous medium. However, it has been generalized to saturated and unsaturated flows in heterogeneous and anisotropic media. q = −K∇h (3) where q denotes specific discharge (volume of fluid per unit cross-sectional area of porous medium per unit time, m/s) and is also called Darcy velocity. K is hydraulic conductivity (m/s), and considered as isotropic, i.e. constant in all directions x, y, z. C. J. Itoua-Tsele, A. F. Ouassissou / Eur. J. Pure Appl. Math, 16 (4) (2023), 2169-2186 2173 Figure 1: Components of hydraulic head in a portion of a groundwater water system. ∇h is the hydraulic gradient, which is the driving force of groundwater flow per unit weight of groundwater (dimensionless). 3.2.1. Isotropic porous medium Darcy’s law for single phase fluid flow in homogeneous isotropic porous media, saturated in terms of total pressure P = p+ ρ0gz (4) with the vertical z axis pointing upward, is written : q = ϕu = −k µ ∇P = −k µ ∇(p+ ρ0gz) = −ρ0gk µ ∇ ( p ρ0g + z ) (5) Where k is the permeabilities (m2), or the intrinsic permeability, a property of the porous medium; ϕ is the porosity of the porous medium (dimensionless), u is the interstitial velocity (m/s), µ is dynamic viscosity (kg/m/s) of the groundwater; p is fluid pressure (Pa), ρ0 is fluid density (kg/m3) and g is the gravitational acceleration (m/s2). Substituting equation (2) into equation (5), we therefore obtain the Darcy’s law in term of hydraulic head h q = KJ = −K∇h (6) where K = ρ0gk µ is hydraulic conductivity of the fluid in a porous medium. Then K is a scalar, and the specific flux vector q with the components qx, qy, qz in the direction of the C. J. Itoua-Tsele, A. F. Ouassissou / Eur. J. Pure Appl. Math, 16 (4) (2023), 2169-2186 2174 Cartesian x, y, z direction and J = −∇h is hydraulic gradient with components Jx = −∂h ∂x , (7) Jy = −∂h ∂y , (8) Jz = −∂h ∂z . (9) Homogeneous isotropic porous medium, K is a constant, so qx = KJx = −K ∂h ∂x , (10) qy = KJy = −K ∂h ∂y , (11) qz = KJz = −K ∂h ∂z . (12) For three dimension flow through inhomogeneous isotropic media we haveK = K(x, y, z). Thus, Darcy’s law in a inhomogeneous domain is: q = −K(x, y, z)∇h (13) In isotropic porous media the off-diagonal elements in the permeability tensor are zero, Kij = 0 for i ̸= j and the diagonal elements are identical, Kii ̸= 0, where subscripts i and j stand for xi, xj , respectively, with x1 ≡ x, x2 ≡ y, and x3 ≡ z. In three dimension space K = Kxx 0 0 0 Kyy 0 0 0 Kzz  (14) In two dimension space K = [ Kxx 0 0 Kyy ] (15) 3.2.2. Anisotropic porous medium In a compressible fluid under isothermal conditions, where, ρ = ρ(p), the pressure head is expressed by ∫ pf p0 dp gρ(p) (16) where p0 is some reference pressure, indicating that the pressure energy stored in the fluid per unit weight of the latter is obtained from the work done in compressing fluid. For such a fluid, it is common to define a piezometric head, h∗, often called Hubbert’s potential. h∗ = h∗(x, t) = z + ∫ pf p0 dp gρ(p) (17) C. J. Itoua-Tsele, A. F. Ouassissou / Eur. J. Pure Appl. Math, 16 (4) (2023), 2169-2186 2175 where z represents the elevation head, that is, potential energy per unit weight of water. The sum of the pressure head and the elevation head is the piezometric head h. Note the use of x in equation (17), representing the coordinates x, y, z to denote a point in a 3D space. Thus, Darcy’s equation is given by the relation q = −k µ [∇p+ ρg∇z] (18) In case of anisotropic porous medium, Darcy’s law is given by q = −K · ∇h (19) where K is tensor of hydraulic conductivity of anisotropic medium. 3.2.3. Principal Directions (i) Although the hydraulic conductivity tensor, K, at a point within an anisotropic porous medium is independent of the coordinate system used, the magnitude of each Kij-component does depend on the chosen coordinate system. Texts on tensor analysis give the rules for transforming these components from one coordinate system to another. (ii) These texts also prove that it is always possible to find three mutually orthogonal directions in space such that when these directions are chosen as the coordinate system for expressing the components, we find that Kij = 0 for i ̸= j (20) Kij ̸= 0 for i = j (21) These directions in space are called principal directions of the hydraulic conductivity tensor (of the anisotropic porous medium). In a heterogeneous porous medium domain, the principal directions may vary from point to point. When the principal directions are aligned with the selected coordinate system, then we have, in three-dimensional space K = Kxx 0 0 0 Kyy 0 0 0 Kzz  (22) and in two dimension space K = [ Kxx 0 0 Kyy ] (23) Thus, in the case of an anisotropic porous medium, if the coordinate axes are oriented along the main axes of anisotropy, the specific discharges along the 3 axes are: C. J. Itoua-Tsele, A. F. Ouassissou / Eur. J. Pure Appl. Math, 16 (4) (2023), 2169-2186 2176 qx = KxxJx = −Kxx ∂h ∂x , (24) qy = KyyJy = −Kyy ∂h ∂y , (25) qz = KzzJz = −Kzz ∂h ∂z . (26) Because the hydraulic conductivity is related to the permeability by the scalar factor K = ρgk µ , the permeability of an anisotropic porous medium is also a second rank tensor. In fact, it is the tensorial nature of the permeability that determines the tensorial nature of the hydraulic conductivity. Sometimes the orientation of the anisotropy directions varies along the hydrogeological structure and thus it is not possible to properly choose the axes of the coordinate systems to satisfy the above-mentioned requirements. The components of the tensor K in a three- dimensional space, can be written in the matrix form K =  Kxx Kxy Kxz Kyx Kyy Kyz Kzx Kzy Kzz  (27) and in a two-dimensional space, as K = [ Kxx Kxy Kyx Kyy ] (28) The hydraulic conductivity tensor is symmetric, that is Kxy = Kyx, (29) Kxz = Kzx, (30) Kyz = Kzy. (31) This means that actually only six distinct components are needed to fully define the hydraulic conductivity in a three-dimensional domain (and only three in a two-dimensional one). Furthermore, the coefficients are non-negative. The specific discharge given by the relation (18) has the following components: qx = −kxx µ ∂p ∂x − kxy µ ∂p ∂x − kxz µ ( ρg + ∂p ∂z ) qy = −kyx µ ∂p ∂x − kyy µ ∂p ∂x − kyz µ ( ρg + ∂p ∂z ) qz = −kzx µ ∂p ∂x − kzy µ ∂p ∂x − kzz µ ( ρg + ∂p ∂z ) (32) C. J. Itoua-Tsele, A. F. Ouassissou / Eur. J. Pure Appl. Math, 16 (4) (2023), 2169-2186 2177 Finally, when the axes of coordinates are oriented along the main axes of anisotropy, the relations (32) become:  qx = −kxx µ ∂p ∂x qy = −kyy µ ∂p ∂y qz = −kzz µ ( ρg + ∂p ∂z ) (33) 3.3. Validity of Darcy’s law Darcy’s law is only applicable when groundwater flow is laminar [8]. The Reynolds number can be used to determine if the flow is laminar or turbulent. RN = ρvd µ (34) where, RN is the Reynolds number, unitless ρ is the fluid density, kg/m3 v is the velocity of fluid, m/s d is the diameter through which moves, m µ is the dynamic viscosity, Pa s. When 1 < RN < 10 flow is considered to be laminar. 3.4. The basic continuity equation In order to derive the continuity equation, which represents the flow of a single phase fluid (water) at constant density in a continuous porous medium under Darcy’s law, [1] used Euler’s approach and assumed that close to the coast, the aquifer is divided into two sub-aquifers, unconfined and confined aquifers; flow always takes place in a three- dimensional in a small control volume that they consider in Cartesian coordinates x, y, z, (see Figure 2). Based on mass conservation law and Darcy’s equation (6), [2, 3] and [1] developed the mathematical equation describing the unsteady state groundwater flow in heterogeneous and anisotropic confined aquifer used in MODFLOW: S0 ∂h ∂t = ∇ · (K∇h) +G, (35) where S0 is the specific storage (m−1), h is the hydraulic head (m), K is the hydraulic conductivity tensor (m/s) and G is the sink/source term (m−1) of water, where negative values are extractions, and positive values are injections and t is time (s) . C. J. Itoua-Tsele, A. F. Ouassissou / Eur. J. Pure Appl. Math, 16 (4) (2023), 2169-2186 2178 Figure 2: Representative elementary volume in a confined aquifer 4. Vertically integrated equations of confined aquifer in two-dimension 4.1. Deriving 2-D Diffusivity Equations by Integration We use the above assumptions to reduce the equation (35) to S0 ∂h ∂t = ∂ ∂x ( Kxx ∂h ∂x ) + ∂ ∂y ( Kyy ∂h ∂y ) +G (36) Now, consider a confined aquifer of variable variable thickness, B (see Figure 3), such that B(x, y) = B2(x, y)−B1(x, y) (37) where B1(x, y) and B2(x, y) are the elevations of its fixed bottom and ceiling, since the aquifer is horizontal and compressible. The Figure 3 shows this aquifer. C. J. Itoua-Tsele, A. F. Ouassissou / Eur. J. Pure Appl. Math, 16 (4) (2023), 2169-2186 2179 Figure 3: Flow through a confined aquifer The vertical integration of the equation (36) along the vertical involving Dupuit’s approximation in (x, y)− plane gives∫ B2(x,y) B1(x,y) S0 ∂h ∂t dz = ∫ B2(x,y) B1(x,y) ∂ ∂x ( Kxx ∂h ∂x ) dz+ ∫ B2(x,y) B1(x,y) ∂ ∂y ( Kyy ∂h ∂y ) dz+ ∫ B2(x,y) B1(x,y) Gdz (38) Now, if the integral is defined as follows∫ b(x) a(x) f(x, t)dt,−∞ < a(x), b(x) < ∞ (39) The derivative of this integral is derived using Leibniz’s integral rule as follows: ∂ ∂x (∫ b(x) a(x) f(x, t) dt ) = ∫ b(x) a(x) ∂ ∂x f(x, t) dt+ f(x, b(x)) · ∂ ∂x b(x)− f(x, a(x)) · ∂ ∂x a(x) (40) Applying this rule on equations (38), we get∫ B2(x,y) B1(x,y) S0 ∂h ∂t dz = ∂ ∂t S0Bh̃, (41) where h̃(x, y) = 1 B ∫ B2(x,y) B1(x,y) h(x, y, z) dz (42) is the average piezometric head along a vertical line at point (x, y).∫ B2(x,y) B1(x,y) ∂ ∂x ( Kxx ∂h ∂x ) dz = ∂ ∂x ( B ˜ Kxx ∂h ∂x ) −Kxx ∂h ∂x ∣∣∣∣ B2 · ∂B2 ∂x +Kxx ∂h ∂x ∣∣∣∣ B1 · ∂B1 ∂x , (43) C. J. Itoua-Tsele, A. F. Ouassissou / Eur. J. Pure Appl. Math, 16 (4) (2023), 2169-2186 2180 and∫ B2(x,y) B1(x,y) ∂ ∂y ( Kxx ∂h ∂y ) dz = ∂ ∂y ( B ˜ Kyy ∂h ∂y ) − Kyy ∂h ∂y ∣∣∣∣ B2 · ∂B2 ∂y + Kyy ∂h ∂y ∣∣∣∣ B1 · ∂B1 ∂y (44) where ˜ Kxx ∂h ∂x (x, y) = 1 B ∫ B2(x,y) B1(x,y) Kxx ∂h ∂x dz (45) and ˜ Kyy ∂h ∂y (x, y) = 1 B ∫ B2(x,y) B1(x,y) Kyy ∂h ∂x dz (46) Thus, the equation (38) gives S0B ∂ ∂t h̃ = ∂ ∂x ( B ˜ Kxx ∂h ∂x ) − Kxx ∂h ∂x ∣∣∣∣ B2 · ∂B2 ∂x + Kxx ∂h ∂x ∣∣∣∣ B1 · ∂B1 ∂x + ∂ ∂y ( B ˜ Kyy ∂h ∂y ) − Kyy ∂h ∂y ∣∣∣∣ B2 · ∂B2 ∂y + Kyy ∂h ∂y ∣∣∣∣ B1 · ∂B1 ∂y +BG̃ (47) According to Dupuit assumption, the hydraulic head is constant along the vertical, then the flow is essentially horizontal , that is, h = h(x, y) and , h (x, y,B1) ≃ h (x, y,B2) ≃ h̃(x, y) leads to ∇h̃ = ∇h, in other word, the gradient of the average head is equal to the average of the head gradient. Therefore, the equation (47) becomes S0B ∂h ∂t = ∂ ∂x [ Txx ∂h ∂x ] + ∂ ∂y [ Tyy ∂h ∂y ] +BG̃ (48) where Txx = K̃xxB and Tyy = K̃yyB (49) are transmissivity components along the x and y directions. Hence, the storage coefficient is S = S0B, the sink/source terme is N = G̃B and we obtain the diffusivity equation in two-dimension for an inhomogeneous isotropic confined aquifer: S ∂h(x, y, t) ∂t = ∇ · [T (x, y)∇h(x, y, t)] +N(x, y, t), (x, y) ∈ Ω ⊂ R2, t > 0. (50) where T (x, y) = ( Txx 0 0 Tyy ) is the aquifer transmissivity (m2/s); h(x, y, t) is the hydraulic head (m), of confined aquifer; S = S(x, y) is the storage coefficient of the C. J. Itoua-Tsele, A. F. Ouassissou / Eur. J. Pure Appl. Math, 16 (4) (2023), 2169-2186 2181 confined aquifer (m−1) and N = N(x, y, t) is the vertical integration of sink/source term, (m/s). The boundary and initial conditions are: Boundary of prescribed head : h(x, y, t) = f1(x, y, t) on Ω (51) Boundary of prescribed flux along such boundary: ∇h(x, y, t) · n(x, y) = f2(x, y, t) (52) where n(x, y) is the outward unit vector normal to the boundary ∂Ω. Initial condition: h(x, y, 0) = f0(x, y) (53) where f0, f1 and f2 are known functions. 5. The Finite Volume Method Now, we will describe the fundamental steps in the derivation of the finite volume method for an unstructured mesh. The principle of finite volume method discretizes the governing equations the computational domain Ω is first split or discretized into ncells smaller control volumes Vi such that the collection of all those subdomains forms a partition of Ω, that is: (i) each Vi is an open, simply connected, these cells may be of arbitrary shape and size, although, traditionally, the cells are convex polygons (in 2D) or polyhedrons (in 3D), i.e., they are bounded by straight edges (in 2D) or planar surfaces (in 3D). (ii) Vi ∩ Vj = ∅ (i ̸= j), (iii) ncells⋃ i=1 Vi = Ω cells are names given to these control volumes, faces are the surfaces that enclose the cells and the vertices of the cells are referred to as nodes, which are linked points within the computational domain. As a result, the governing equations are integrated over each control volume Vi, and the resultant integrals are approximated with the mid-point rule or Gaussian quadrature rule [6, 12]. The volume-averaged value is defined as hO = 1 VO ∫ V0 h dV (54) Theorem 1 (Gauss’s divergence theorem). Let us consider a closed bounded region in space of volume VO bounded by a surface of area S, as shown schematically in Figure 4. C. J. Itoua-Tsele, A. F. Ouassissou / Eur. J. Pure Appl. Math, 16 (4) (2023), 2169-2186 2182 Let J be a vector field that is continuous and has continuous first partial derivatives in this region. Gauss’ divergence theorem asserts that∫ VO ∇ · J dV = ∫ S J · n̂ dA (55) where n̂ is the outward-pointing unit surface normal to the differential area dA. Figure 4: Control Volume [12] In general, finite volume approaches may be identified by the following criteria [9]: (i) the geometric shape of the control volumes Vi, (ii) the position of the unknowns (“problem variables”) with respect to the control vol- umes, (iii) the approximation of the boundary (line (d = 2) or surface (d = 3)) integrals. The second criterion, in particular, splits the finite volume method into two broad cate- gories: the cell-centered and the cell-vertex finite volume methods. In the cell-centered methods, the unknowns are associated with the control volumes (for example, any control volume corresponds to a function value at some interior point (e.g., at the barycentre)). In the cell-vertex methods, the unknowns are located at the vertices of the control volumes. 6. Discretization of the diffusivity equation To discretize equation (50), we integrate it into each cell VO, and for each time interval [t, t+∆t], yiels:∫ t+∆t t ∫ VO S ∂h ∂t dV dt = ∫ t+∆t t ∫ VO ∇ · (T · ∇h) dV dt+ ∫ t+∆t t ∫ VO N dV dt (56) C. J. Itoua-Tsele, A. F. Ouassissou / Eur. J. Pure Appl. Math, 16 (4) (2023), 2169-2186 2183 6.1. Discretization of Transient term For the transient term, we integrate with respect to t and apply the average definition, then this gives: ∫ t+∆t t ∫ VO S ∂h ∂t dV dt = ∫ VO S(hn+1 − hn) dV = SOVO(hO − holdO ) (57) 6.2. Discretization of Diffusion term ∇ · (T · ∇h) The discretization for the diffusion term is given by∫ t+∆t t ∫ VO ∇ · (T∇h) dV dt = ∫ t+∆t t ∫ S T (∇h · n̂) dAdt = ∫ t+∆t t Nf,O∑ f=1 ∫ Sf T (∇h · n̂) dA  dt = ∫ t+∆t t Nf,O∑ f=1 Tf (∇h · n̂)fAf  dt = Nf,O∑ f=1 Tf (∇h · n̂)fAf∆t (58) The estimation for the flux, (∇h · n̂)f , is given as: (∇h)f · n̂f = hN − hO δf − [ ha(f) − hb(f) δfAf ] (xa − xb)(xN − xO) + (ya − yb)(yN − yO) Af . (59) where ha(f) and hb(f) are the value of h at the vertices a and b, respectively. To compute ha(f) and hb(f), we make use of the interpolation function facteur, if a vertex (or node) is influenced by Nc cells , then wv,i = 1/di Nc∑ i=1 1/di , (60) where the interpolation function facteur wv,i, represents the contribution of the ith sur- rounding cell to the vertex, v and di is the distance between the vertex v and the cell center i of cell i. Once the interpolation function has been computed and stored, the vertex (or nodal) value of h can be computed using hv = Nc∑ i=1 wv,ihi (61) C. J. Itoua-Tsele, A. F. Ouassissou / Eur. J. Pure Appl. Math, 16 (4) (2023), 2169-2186 2184 We chose the inverse weighted function because the cell effect will decrease as the distance increases. 6.3. Discretization of Source Term In this study, we consider that the sources take the form of point sources and sinks N i(x, t) which is located at points xi and at time t, we express the source term in the form [3]: N = Ns∑ i=1 N i (xi, t) δ (x− xi) . (62) where Ns is the number of sources and sinks and δ(x−xi) denotes the Dirac delta-function for two dimensions, defined formally by δ (x− xi) = lim a→0 { 1/a2 |x− xi| < a 0 elsewhere (63) Actually, N i can represent both pumping wells and recharging (with N i < 0 and N i > 0 respectively). The discretization for the source term is given by,∫ t+∆t t ∫ VO N dV dt = ∫ t+∆t t ∫ VO Ns∑ i=1 N i (xi, t) δ (x− xi) dV dt = Ns∑ i=1 N i (xi, t)O δ (xO − xi)VO∆t (64) where the subscript O indicates the cell center of the control volume O. 6.4. Numerical Schemes In this study, for stability reasons, we use the implicit approach. Combining Equations (57), (58) and (64), the complete discretization for equation (50) can now be written for each VO as SOVO(hO − holdO ) = Nf,O∑ f=1 Tf ( hN − hO δf − [ (∇h)f · t̂f ] δf t̂f · 1f ) Af∆t + Ns∑ i=1 N i (xi, t)O δ (xO − xi)VO∆t (65) Finally, equation (65) is assembled in matrix form in order to obtain the discrete equations between each control volume VO value and its neighboring values as follows: aOhO = Nf,N∑ f=1 af,Nhf(N) + bN (66) C. J. Itoua-Tsele, A. F. Ouassissou / Eur. J. Pure Appl. Math, 16 (4) (2023), 2169-2186 2185 where the coefficients are given by aO = SOVO + Nf,O∑ f=1 Tf hO δf Af∆t, (67) af,N = Tf δf Af∆t, (68) bN = Ns∑ i=1 N i (xi, t)O δ (xO − xi)VO∆t− Nf,O∑ f=1 Tf [ (∇h)f · t̂f ] δf t̂f ·1fAf∆t+SOVOh old O . (69) 7. Conclusion The main objective of this study was to develop the vertically integrated groundwater flow equation, under Dupuit’s assumption, the diffusivity equation was then derived, the study area is often complex, and it required the relevant numerical method, finite volume method has the advantage of working on complex geometric and it is also stable, the implicit FVM schemes were developed for 2D. The obtained numerical scheme will be used to evaluate the distribution of the hydraulic head in the study area and investigate the fluctuation of groundwater level, from which the drawdown phenomenon can be observed. 8. Recommendations for further research Based on the findings of this study, the researcher found that there are still some areas of mutual interest. The researcher now suggests the following recommendations for future studies as a result of this work. • Unsteady state solution of the confined anisotropic and heterogeneous aquifers to two-dimensional flow • Validate the developed models by comparing them with other methods Conflict of interest The authors of this article acknowledges its an original work from their investigation and therefore there is no conflicts of interest whosoever as regards its publication. Acknowledgements The authors wish to express their gratitude to Prof. Christian TATHY and Prof. David Mumo Malonza for their contributions to this study who so kindly donated their time and insights. REFERENCES 2186 References [1] Mary P Anderson, William W Woessner, and Randall J Hunt. Applied groundwater modeling: simulation of flow and advective transport. Academic press, 2015. [2] Jacob Bear. Hydraulics of groundwater. Courier Corporation, 2012. [3] Jacob Bear and Alexander H-D Cheng. Modeling groundwater flow and contaminant transport, volume 23. Springer, 2010. [4] Prosun Bhattacharya, AL Ramanathan, Arun B Mukherjee, Jochen Bundschuh, D Chandrasekharam, and AK Keshari. Groundwater for sustainable development: problems, perspectives and challenges. 2008. [5] Huseyin Y Dalkilic and Amin Gharehbaghi. Numerical modeling of groundwater flow based on explicit and fully implicit schemes of finite volume method. 2021. [6] Robert Eymard, Thierry Gallouët, and Raphaèle Herbin. Finite volume methods. Handbook of numerical analysis, 7:713–1018, 2000. [7] Bahareh Kalantar. Groundwater Management and Resources. BoD–Books on De- mand, 2021. [8] Michael Kasenow. Applied ground-water hydrology and well hydraulics. Water Re- sources Publication, 2001. [9] Peter Knabner and Lutz Angermann. Numerical methods for elliptic and parabolic partial differential equations. 2021. [10] Makosha Ishmaeline Charlotte Morakaladi and Abdon Atangana. Mathematical model for conversion of groundwater flow from confined to unconfined aquifers with power law processes. Open Geosciences, 15(1):20220446, 2023. [11] Rachel Rotz. Hydrogeologic properties of earth materials and principles of ground- water flow, 2021. [12] Mazumder Sandip. Numerical methods for partial differential equations: finite differ- ence and finite volume methods. Academic Press, 2015. [13] V Vázquez-Báez, A Rubio-Arellano, Dolores Garćıa-Toral, I Rodŕıguez Mora, et al. Modeling an aquifer: Numerical solution to the groundwater flow equation. Mathe- matical Problems in Engineering, 2019, 2019. [14] Ming-Chang Wu and Ping-Cheng Hsieh. Variation of groundwater flow caused by any spatiotemporally varied recharge. Water, 12(1):287, 2020. [15] Xiaoping Zhang, Yan Zhu, Jing Wang, Lili Ju, Yingzhi Qian, Ming Ye, and Jinzhong Yang. Gw-pinn: A deep learning algorithm for solving groundwater flow equations. Advances in Water Resources, 165:104243, 2022.