# Numerical Simulation of a Two-Dimensional Groundwater Pollute Transport Problem Using Incompressible Steady-State Navier-Stokes Equations and Diffusion-Convection Equations

Authors: Jacob Nyende; Isaac Enyogoi; John Mango; Henry Kasumba
Journal: Modelling and Simulation in Engineering
 (2022)

---

## Abstract

Most of the real contaminant problems are defined domains that are geometrically complex and can have different boundary conditions in different areas. Therefore, it is usually difficult to find a solution analytically, so we use the approximate method to generate an approximate function. One answer to this problem is the finite element approach (FEM). This study presents a partial differential equation (PDE) simulation system that uses numerical techniques for the distribution of pollutant concentrations in groundwater in space and time. The movement of the liquid is described by the incompressible steady-state Navier-Strokes equation, while the transport of pollutants is described by the diffusion-convention equation. The variation formulation that forms the basis of FEM and MATLAB is discussed along with the selection of the abstract approximation space and the welfare of the weak formulation. The motivation for this study comes from a specific and considered water body with the discharge of factory effluents on the ground that ends up reducing the quality of groundwater. First, the fluid flow equation is solved to obtain velocity and pressure profiles. Steady-state concentration profiles were obtained for various values of diffusion coefficient (D), baseline, and input concentrations. The results showed that decreasing the diffusion coefficient D increased the number of pollutants for convective transport and decreased the number of pollutants that diffused from the entrance. Although groundwater is not completely safe, it is concluded that experimental studies are necessary decision-making basis for water resource protection, especially in water pollution emergencies.


---

## Body

## 1. Introduction

In this paper, a model based on a stationary incompressible 2-dimensional Navier-Stokes equation for fluid flow and a 2-dimensional convection-diffusion equation for polluting water transport (one type of pollutant) is then formulated and the numerical model solved for different boundary conditions and, finally, simulations to predict the movement and concentration of the pollutant in regions occupied via of porous media at different positions and times. The driving motivation for this work comes from a specific and considered industrial problem in the transport and dispersion processes of polluted particles that end up reducing the quality of groundwater [1]. In the literature, water is a basic need of life, used domestically in industry, recreation, and irrigation purposes. Groundwater is found in underground aquifers and provides approximately 97% of the world’s consumable water [2]. In Uganda, many communities, especially the rural poor, depend on streams and swamps for supply of water [3]. In recent years, the emphasis for policymakers and researchers has shifted from water accessibility to water quality problems. Environmental pollution problems require prompt action to prevent the reduction of water quality. The issue of water quality is an important factor in sustainable human development [4]. Groundwater is polluted due to anthropogenic activities on the land, which include the use of agrochemicals and poor disposal of domestic and industrial waste. Increasing trends in the use of fertilizers and pesticides in agricultural production systems and disposal of industrial waste are the reasons for contamination of groundwater [4]. A study on risk factors contributing to microbiological contamination of near subsurface water in Kampala indicated that water from shallow protected springs was polluted [5]. Pollution level increased in the rainy season due to storm water runoff and infiltrated into the groundwater [6]. The impact of pollution on shallow groundwater in Bwaise III Parish, Kampala [7], showed contamination of groundwater; in addition, the pollution levels of water resources in Pece Wetland, Gulu Municipality, indicated that domestic water sources were contaminated [8]. Modelling fluid flow problems falls in the branch of applied mathematics called computational fluid dynamics (CFD), which uses tools from numerical analysis, fluid dynamics, and computer science. Discretization techniques which include boundary element, finite difference, finite volume, and finite element methods have been globally utilized in solving solute transport and distribution of fluids [9]. Problems solved by these discretization techniques include multiphase and single-phase fluid flows in porous medium [10–12], flow in a thermal labyrinth [13], pollute transport in the atmosphere [14, 15], pollute transport in rivers [16, 17], and pollute the groundwater [2]. FEM is a numerical analysis technique which employs the concept of piecewise “approximation” to approximate partial solutions of differential equations. It is a powerful discretization tool which can accurately discretize a domain of any size or shape. The Navier-Stokes equations are fundamental models in fluid dynamics [18] that describe motions of fluids in ℝn,n=2or3. These equations solve velocity vector uix,t1≤i≤n∈ℝn and pressure p∈ℝ for any position x∈ℝn and time t≥0, [19]. Many groundwater problems have been modelled by Darcy’s law [1, 12, 20–22], among others. Henry Darcy discovered that water flow through a porous medium was governed by equation (1).
(1)Q=KAμ∇p−ρgz∇l,whereQ is the total discharge through the medium, K is the permeability of the medium, A is the cross-sectional area of the flow, ∇p−ρgz is the pressure gradient, μ is the viscosity, and ∇l is the length over which pressure drop occurs [21].In this study, fluid flow is modelled by the Navier-Stokes equations instead of the usual Darcy’s law because, under certain assumptions, Darcy’s law is derived from the Navier-Stokes equations, [23]. The momentum balance (equation (2)) can be written as
(2)ρDvDt=∇.τ+ρF+T,where by:Dv/Dt=∂v/∂t+v.∇v is the material time derivative of a given particle, τ is the stress tensor, F is the external body force per unit mass, and T is the drag force/unit volume exerted on the fluid. The momentum balance (equation (1)) for the flow of water in the saturated zone reduces to the well-known Darcy’s law under certain conditions [23]. Moreover, Darcy’s law is reliable for values of Reynolds’number<1. At very low Reynolds’ number, the circulation generated can be neglected, but as Reynolds’ number increases, it becomes important and a noticeable contribution to the resistance to flow. The fact is that the circulation zone becomes smaller at even higher Reynolds’ number. Pollutants generally dissolve and are carried by infiltrating rain water into unsaturated soil above the water table. Pollutants then enter the saturated zone and migrate in the direction of groundwater flow (Figure 1).Figure 1
Schematic diagram illustrating groundwater contamination from a waste site (source: Environment Canada, researchgate.net, accessed November 4, 2021).


## 2. Materials and Methods

Human health is harmed by suspended particles in water. Particle transportation and distribution of matter suspended in water is associated with fluid motion and turbulence. CFD is the most appropriate modelling approach for studying particle spatial distribution in a particular domain [9]. These methods enable complex and abstract mathematical models or theories to be easily visualized through the use of simulations. In the formulation of groundwater models, errors can arise from conceptual problems due to excluding relevant or including irrelevant physical process in the model or using an inapplicable model for the groundwater solute transport problem and numerical errors like truncation and round-off errors and using wrong input data [1]. The CFD model is determined by the nature of the physical process to be simulated, the objective of the study, and the available resources [24]. Equation parameters can easily be adjusted to observe how these changes will affect the model [25]. Because the mathematical model is the result of coupling multiple physical fields, which restricts the choice of function spaces when using the finite element method, an adequate choice of function spaces for pressure and velocity is required [26]. Mesh generation depends on the complexity of the geometry of the domain and can either be structured, block-structured, or unstructured. For example, in a one-dimensional case, if the computational domain is divided into N equal subintervals, we generate a structured mesh, in the case of a complex domain, that requires a fully unstructured mesh. A good and reliable numerical model should be able to simulate a transport phenomenon accurately and suppress the instability arising out of discretization in the computational domain [27]. Groundwater flow and solute transport are described by second-order partial differential equations that can be parabolic, elliptic, or hyperbolic. This classification is essential because the choice of a numerical method should be best suited to the nature of the PDE. For instance, a diffusion-convection process dominated by the convective term approximates a hyperbolic type of equation, while a diffusion-convection process dominated by the diffuse term approximates a parabolic type of equation. FEM is a universal discretization tool for partial differential equations. Its advantage is that it can easily be defined on a complicated geometry and can improve the quality of the numerical solution by fine-tuning the model grid. The mesh need not be uniform; it can be made finer in areas of the domain where greater accuracy is needed. Its downside, however, is that the programming of FEM is more complicated and needs standard software packages. The finite difference method (FDM) is simpler and easier to program, but it is only applicable to a regular domain [28]. The extract from Konikow and Reilly [28] is an illustration of the finite difference and finite element methods in discretizing a domain of a field. The rectangular grid of the FDM fails to cover the entire domain of the field, while the triangular grid of the FEM can easily be constructed as shown in Figure 2.Figure 2
Hypothetical application of finite difference (b) and finite element (c) grids to an irregularly bounded aquifer (a) (adopted and modified from Konikow and Reilly [28]).
(a)(b)(c)When using the FDM technique, the domain is divided into uniform subareas and the PDE’s are then replaced by approximating functions which are written in terms of finite subspaces. The time derivative is also approximated by a finite difference expression obtained using a Taylor series expansion [23]. When formulating the finite element method, the variational direct energy balance or weighted residual approaches can be used. The weighted residual and variational approaches are the most commonly used approaches for groundwater models [1]. When using the variational approach under FEM, the appropriate PDE together with its boundary conditions and initial conditions are replaced with a corresponding variational problem. In this approach, the unknown function uh is approximated throughout the solution domain as in Equation (3)
(3)uhx,y=∑i=1n−1ciϕjx,y,whereϕj are suitable bases or shape functions, coefficients ci are unknowns, and n is the number of nodes.A finite element Newton method for the solution of steady-state Navier-Stokes equations for 2-dimensional incompressible flows was discussed [29]. In this work, a weak variational formulation of the problem formulation and an unequal order interpolation for pressure and velocity were adopted. A Newton method was used for the nonlinear system of coupled equations written in an incremental form and the Jacobian linear system was solved using a direct algorithm. They explained that the Newton iterative method was a suitable technique because of its efficiency since only a few iterations were sufficient for convergence to a very accurate steady solution provided the initial guess was not chosen too far from the solution. They also explained the use of the incremental form of the Newton iteration, which fully exploited the quadratic nature of the nonlinearity of the Navier-Stokes system and thereafter obtained an algorithm which was optimal both in the imposition of the boundary conditions and with respect to the cancellation of numerical errors. For solutions at high Reynolds’ numbers, the Newton method may fail when the Stokes flow is too far from the nonlinear solution. They then resorted to a continuation stepwise manner by computing different intermediate solutions for smaller values of the Reynolds’ numbers. The Laplace transformation technique has been commonly used to obtain the solution of the diffusion-convection equation because of its simplicity compared to other methods [30]. Other analytical solutions to the convection-dispersion solute transport equation are discussed [31]. A groundwater contamination, diffusion, and spreading model can be solved by Fourier transforms [2]. A similar model could have been done without assuming that all covering parameters are constant because this is highly unlikely in reality. For example, soil conductivity, which measures how fast a fluid can be transported through the soil, is a highly variable quantity. Other models of pollution in river type aquatic systems can be found [4]. To prevent groundwater contamination in both shallow and deep aquifers [32], cut-off walls are often used. Discacciati and Quarteroni [12] utilized a model that couples the Navier-Stokes and Darcy equations and modelled the filtration of incompressible fluids through a porous medium [33]. In this work, the motion of incompressible fluids was described by the Navier-Stokes equations while the filtration process was described by Darcy’s equation. A computational domain divided into two parts was considered, one for the fluid and the other for the porous media, and went ahead to discuss the coupling conditions (Beavers-Joseph-Saffman) including their mathematical justification via a homogenization theory. Interface conditions were also introduced because of the multidomain structure of the model. In this regard, a much better model for pollute transport in a groundwater contaminant problem would require the coupling of the Navier-Stokes equations with a groundwater flow equation in the porous media together with a diffusion-convection equation for pollute transport in the two regions.


### 2.1. A PDE Model for Groundwater Flow

To achieve the objectives of groundwater flow and pollute transport problem, a groundwater flow equation based on an incompressible steady 2-dimensional Navier-Stokes equation and a 2-dimensional convection-diffusion equation describing pollute transport is constructed. The Navier-Stokes equations are derived on the basis of the conservation principles of mass (continuity) and linear momentum (Newton’s second law of motion). The first law of thermodynamics (energy) is not considered because of the assumption that heat transfer in groundwater is negligible. Finally, the boundary conditions both Dirichlet and Neumann are stated, together with reasonable assumptions. The density of water is assumed constant, that is, it is not affected by variations in concentration. The pollute transport equation is derived based on a Fickian model which assumes that the waste is transported mainly due to the concentration gradient and that the dispersive flux occurs in a direction from higher concentrations towards lower concentrations [34]. In this work, the PDE model for pollute transport in groundwater was formulated by first deriving an equation for groundwater flow using Navier-Stokes equations.


### 2.2. Derivation of the Navier-Stokes Equation

The flow of water is described using steady incompressible Navier-Stokes equations. These equations are derived on the basis of the conservation principles of mass (continuity) and linear momentum. LetΩ∈ℝ2 describe the domain covered by the flowing water. A fluid particle X (Figure 3) is considered to be moving through x∈ℝ2 at time t.Figure 3
Mapping from the Langrangian coordinate to Eulerian coordinate.The Lagrangian or material coordinate isX. It is the observed trajectory of a specific fluid particle as the fluid flows while x is a Eulerian coordinate because the fluid particle is observed through a fixed region of interest [35].Also, letωo⊂Ω be a subset of t=0. Define the function ϕ:ℝ2⟶ℝ2 such that det ∇xϕX>0 for all X∈Ω the change of the particle’s position is described by equation (4).
(4)ωt=ϕX,t:X∈ω0=ϕω0,t.Defineu=ux,t to be a vector field on Ω at time t, ρx,t the mass density, ρRX,t the mass density of the reference frame f¯x,t the force density, tx,t,n the surface force density also known as traction or contact forces, and n the outward vector to Ωt. By the conservation of mass property, the total mass remains constant as shown in equation (5).
(5)∫ω0ρRXdx=∫ωtρx,tdx=constantforallt>0.Since the total mass remains constant, the rate of change of mass inωt is given by
(6)ddt∫ωtρx,tdx=0.To simplify the term on the left-hand side of equation (6), Reynolds’ transport (Theorem 1) leading to equation (7) was used.Theorem 1.
For an arbitrary single-valued scalar functionρ=ρx,t with continuous derivates and an arbitrary control region ωt with boundary ∂ωt, outward-pointing unit-normal n, and u the local velocity of the fluid across ∂ωt, the following integral relation holds:
(7)ddt∫ωtρx,tdx=∫ωt∂ρx,t∂tdx+∫∂wtρx,tu.nds,
whereds is an element of length on ∂ωt.Source: Reinstra and Hirschberg [36] and Powers [37].In simple terms, the rate of change of a quantity,ρ=ρx,t, in ωt is a result of the rate at which the quantity is being produced in ωt plus the net amount that crosses the boundary ∂ωt Reynolds’ transport theorem is used to translate the integral conservation laws into differential conservation laws. This helps in formulation of the basic conservation laws of fluid dynamics. To simplify calculations, the boundary integral on the right-hand side of equation (7) is transformed to the domain integral using Gauss divergence (Theorem 2).Theorem 2.
Letsωt a 2-dimensional domain with boundary ∂ωtds an element of length (arc length) on ∂ωt and associated unit outward normal vector n with u a vector field on ωt then ∫ωt∇.udx=∫∂wtu.nds, source: [36]. Gauss’ divergence theorem 1.1 shows the relationship between a domain and surface integral. Applying the Gauss divergence theorem, the right-hand side of equation (7) can be written as
(8)∫∂wtρx,tu.nds=∫ωt∇.ρx,tudx.
Substituting equations (8) into (7) yields equation (9) as
(9)ddt∫ωtρx,tdx=∫ωt∂ρx,t∂tdx+∫ωt∇.ρx,tudx,=∫ωt∂ρx,t∂t+∇.ρx,tudx.
Applying equation (6) to equation (9),
(10)⇒∫ωt∂ρx,t∂t+∇.ρx,tudx=0.Since the integral over ωt is arbitrary, equation (10) can be written as
(11)∂ρx,t∂t+∇.ρx,tu=0Equation (11) is known as the continuity equation. Applying the incompressibility condition ρx,t=ρ=constant to equation (11), the continuity equation now becomes ∇.u=0. The conservation of momentum equation was derived from Newton’s second law of motion. From f¯=ma=d/dt∫ωtρx,tdxu=d/dt∫ωtρx,tudx, where f¯ is the force action on the particle, m is the mass, a is the acceleration, and u is the velocity of the particle. The total force f¯ acting on the subdomain ωt is the sum of the field force fx,t and the surface force tx,t,n.By Newton’s second law of motion,d/dt∫ωtρudx=∫ωtρfx,tdx+∫∂wttx,t,n.nds, where ds is an element of length on ∂ωt. Considering the action of the body or field forces fx,t and the normal component of the surface forces tx,t,n=σx,t.n, where σx,t∈ℝ2×2 is the fluid stress tensor, the formula is written as
(12)ddt∫ωtρudx=∫ωtρfx,tdx+∫∂wtσx,t.nds.Considering theith component of the vectors ux,t,f¯x,t, and σi,∗x,t, where σi,∗x,t denotes the ith row of the matrix σx,t, equation (12) is rewritten as
(13)ddt∫ωtρuidx=∫ωtρfidx+∫∂wtσi,∗.nds,fori=1,2.Using the divergence theorem, the boundary integral on the right-hand side of equation (13) is written as
(14)∫∂wtσi,∗.nds=∫∂wt∇.σi,∗dx.Moreover, by Theorem1, the term on the left-hand side of equation (13) can be written as
(15)ddt∫ωtρuidx=∫wt∂ρui∂t+∇.ρuiudx,=∫ωt∂ρui∂t+ui.∇.ρu+∇ui.ρudx.And by the incompressibility condition, the following is obtained:(16)ddt∫ωtρuidx=∫wt∂ρui∂t+u.∇ρuidx, substituting equations (14) and (16) into equation (13),
(17)∫wt∂ρui∂t+u.∇ρuidx=∫ωtρfidx+∫∂wt∇.σi,∗.dx,⟹∫wt∂ρui∂t+u.∇ρui−ρfi−∇.σi,∗dx=0.Since the integral overωt is arbitrary, equation (17) becomes
(18)∂ρui∂t+u.∇ρui−ρfi−∇.σi,∗=0,where∇.σ=∂σ11∂x+∂σ12∂y∂σ21∂x+∂σ22∂y,∇u=∂u1∂x+∂u1∂y∂u2∂x+∂u2∂y.Introducing constitutive relations into equation (18), σ=−pI+τ and p∈ℝ,I∈ℝ2×2; τ∈ℝ2×2=λ∇.uI+2μD, where λ is the volume viscosity, μ is coefficient of shear viscosity, and we set the deformation tensor D in
(19)D=12∇u+∇uT, is the small strain tensor [38]. Substituting for p, D and σ for ∇.σ in equation (18),
(20)∇.σ=∇.−pI+λ∇.uI+μ∇u+∇uT=−∇p+λ+μ∇∇.u+μ∇.∇u,where∇.∇u=∇2u=Δu is the Laplacian of u. Applying the incompressibility condition to equation (20),
(21)⟹λ+μ∇∇.u=0.Therefore, equation (21) is obtained as
(22)∇.σ=−∇p+μ∇u. Substituting equation (22) into equation (18),
(23)⟹∂ρu∂t+u.∇ρu−ρf+∇p−μΔu=0.Assuming the mass densityρ to be constant, equation (23) can be written as
(24)∂u∂t+u.∇u=f−1ρ∇p+vΔu,wherev=μ/ρ is the kinetic viscosity.To make an analysis of the physical problem when the equation parameters of equation (24) change, scaling of the Navier-Stokes equation (24) was carried out. Normalization of these scales resulted in the formulation of dimensionless groups such as Reynolds’, Froude’s, Euler’s, Weber’s, Prandtl’s, and Mach’s numbers which represent the relative importance of various parts of the Navier-Stokes equations and are also key in determining the flow regimes. Scaling also helped to reduce the number of variables. For instance, the combined effect of both viscosity and density can be determined by one dimensionless variable called Reynolds’ number.


### 2.3. Nondimensionalization of Navier-Stokes Equations

Let(25)u~=uU,t~=tT,x~=xL,T=LU,whereU, L, and T are characteristics speed, length, and time, respectively. From expressions in (25), we obtain
(26)∂u∂t=∂Uu~∂Tt~=UT∂u~∂t~=UTut.~In equation (26), the convective termu.∇ubecomes
(27)u.∇u=Uu~.∇Uu~=U2u~.∇u~=U2u~.∂u~∂x~L=U2u~L.∂u~∂x=U2Lu~.∇~u~.The pressure term becomes(28)1ρ∇p=1ρ∂p∂x=1ρ∂p∂x~L=1ρL∂p∂x~=1ρL∂p∂x~=1ρL∇~p.The viscous termvΔu in (29) becomes
(29)vΔu=v∂2u∂x2=v∂2Uu~∂x~L2=vL2∂2Uu~∂x~2=vUL2∂2u~∂x~2=vUL2Δ~u~,(30)∇.u=∇.u~U=U∇.u,~(31)⟹∇.u=U∂∂x.u~=U∂∂x~L.u~=UL∂∂x~.u~,(32)⟹∇.u=UL∇~.u~=0. Substituting equations (26)–(30) into equation (24), we obtain
(33)⟹UTut~+U2Tu~.∇~u~=f−1ρL∇~p+vUL2Δ~u~.Also letp~=1/ρU2P, then
(34)∇~p=∇~ρU2p~=ρU2∇~p.Introducing a nondimensional source term(35)f~=LU2f,and substituting equations (34) and (35) into equation (33) and further substituting T=L/U and multiplying L/U2, we obtain
(36)u~t+u~.∇~u~=f~−∇~p~+1ReΔ~u~.whereRe=LU/v=ρLU/μ is the Reynolds number. It is a dimensionless ratio of inertial to viscous forces used to determine the flow regime. Dropping the tilde notation and assuming steady state, equations (36) and (28) form
(37)−αΔu+u.∇u+∇p=p.And∇.u=0, where α=1/Re, equation (36) represents the dimensionless steady incompressible Navier-Stokes equations.Pollute transport is due to the combined effect of diffusion and advection in a fluid. The PDE describing it as derived from the principle of conservation of mass using Fick’s law [30]. Underground water flows through the pores (space between soil particles) in the soil. Some soils are more compact than others. The ease with which water flows through the soil depends on the soil porosity denoted by ϵ∈0,1 such that
(38)ϵ=TotalvolumeofvoidsinthesoilTotalvolumeofsoil.


### 2.4. Diffusion-Convection Equation

Consider a porous domainΩ∈ℝ2 and ωt⊂Ωt with Cx,t denoting the density or concentration of the pollute in mass per unit area, Qx,t denotes the flux, that is, mass per unit length per unit time crossing the boundary qx,t the rate at which the pollutant is increasing or reducing in ωt due to the source term. The total mass of pollute in ωt is ∫ωtϵCdx and the amount of pollute that crosses the boundary ∂ωt is ∫ωtQ.nds., where ds is an element of length on ∂ωt and n is a unit outer normal vector to ds. The pollutant crosses the boundary ∂ωt in three ways:
(1)
Advection: this is due to the bulk flow of water molecules where pollute particles are carried along the streamlines. It is given byQα=Cv (concentration X velocity)(2)
Molecular diffusion: this is caused by the random motion and collusions of pollutant particles. This occurs even when the fluid mass is at rest(3)
Mechanical dispersion: this is due to the velocity variations caused by the different flow paths that the groundwater takesThe total increase or decrease of the pollutant inωt due to the source terms is given by ∫ωtqdx. By the mass balance principle, the rate of change of the total mass of the pollutant in ωt equals the net rate of pollutant mass that flows through the boundary ∂ωt plus the net rate of increase or decrease of pollutant mass in ωt; this is summarized as
(39)localchange=advection+diffusion+source.Thus,(40)ddt∫ωtϵCdx=−∫ωtQ.nds+∫ωtqdx.By the divergence theorem, the first term on the right-hand side of equation (40) can be written
(41)−∫ωtQ.nds=−∫ωt∇.Qdx,so that equation (41) becomes
(42)ddt∫ωtϵCdx=−∫ωt∇.Qdx+∫ωtqdx.Applying the Reynolds transport theorem, Gauss divergence theorem, and the incompressibility condition, the derivative term on the left-hand side of equation (42) is written as
(43)ddt∫ωtϵCdx=∫ωt∂∂tϵCdx+∫ωt∇.ϵCvdx,which further simplifies to(44)∫ωt∂∂tϵC+ϵv.∇C+∇.Q−qdx=0.Sinceωt is arbitrary, equation (44) implies
(45)ϵ∂∂tϵC+ϵv.∇C+∇.Q−q=0.The contribution to flux from the effect of diffusion (molecular) and solubility (mechanical dispersion) according to Fick’s law:(46)Qd=−Dx,t∇ϵC,whereDx,t is the diffusion coefficient. The equation means that the flux is proportional to the negative gradient of concentration and the pollution will diffuse from a region of high to low concentration. The total flux is the sum of advective, molecular, and mechanical dispersions [34]. Therefore,
(47)Q=Qα+Qd=−Dx,t∇ϵC+Cv,wherevx is the convection velocity also known as advection velocity, substituting equation (47) into equation (45) and assuming the diffusivity Dx,t is constant, ϵ∂/∂tC+ϵv.∇C=−∇.−D∇ϵC+Cv+q, and by the incompressibility condition ∇.v=0, we obtain
(48)ϵ∂∂tC+ϵv.∇C=ϵDΔC−v.∇C+q,which simplifies to(49)∂C∂t=DΔC−u.∇C+F,whereu=ϵ+1ϵ+1/ϵv is the average velocity vector of groundwater obtained from solving equation (37) and F=q/ϵ is the source term. Combining equations (37) and (49), the complete model equations which are a set of highly nonlinear partial differential is as follows;
(50)−αΔu+u∇u+∇p=f,inΩ,∇.u=0,inΩ,∂C∂t−D∇C+u.∇C=F,inΩT≔0,T×Ω.With initial conditionCx,0=C0 and boundary conditions prescribed by
(51)u=gDon∂ΩD−pn+α∇u.n=gNon∂ΩΝ,where∂ΩD and ∂ΩΝ represent the domain boundary prescribed by Dirichlet and Neumann conditions, respectively, and ∂Ω=∂ΩD∪∂ΩΝ.Also, it is expressed as(52)C=CDon∂ΩDand∂C∂n=CNon∂ΩN,whereC0,CD,CN,gD, and gN are known data to be prescribed later. Dirichletian boundary conditions are called essential because they are imposed explicitly on the solution, that is, on the test space, while Neumann boundary conditions are called natural because they are implicitly defined in the variationally weak formulations in which the model is based on the following assumptions:
(a)
It is known a priority regions where groundwater flows(b)
Groundwater flows at a rate which depends on the permeability of the soil(c)
No chemical reactions occur during the interaction of pollutants with waterThe existence of solutions to equation (50) is discussed below.


### 2.5. Existence of Solutions to the Model Equations

The following definitions are used to describe the functional spaces.Definition 3.
The space of real valued functions that are square integral onΩ with respect to the Lebesgue measure dx is denoted by the relationship
(53)L2Ω≔u:Ω⟶ℝ∣∫Ωu2xdx<∞,(54)L2Ω2≔u:Ω⟶ℝ2∣uiϵL2,i=1,2,
The spaceL2Ω2 is equipped with the norm
(55)uL2Ω2=∫Ωux2dx1/2,and scalar products in [39]
(56)u,vL2Ω2=∫Ωux.vxdx.Definition 4.
The Sobolev spaceH1Ω is denoted by the expressions in [39]
(57)H1Ω≔υ:Ω⟶ℝ∣υϵL2Ω,∇υϵL2Ω2,(58)H1Ω2≔v:Ω⟶ℝ2∣uiϵH2,i=1,2.
This space is equipped with the scalar product in(59)u,vH1Ω2=∫Ωuv+∇vdx,that induces the norm in [39]
(60)uL2Ω2=∫Ωu2+∇u2dx1/2.


### 2.6. Functional Spaces

LetV0,ψ,andW0 be test spaces for velocity, pressure, and pollute concentration, respectively, as described below;
(61)V0≔uϵH1Ω2∣u=0on∂ΩD,ψ≔ψϵL2Ωand,W0≔ωϵH1Ω∣ω=0on∂ΩDThe velocity solution spaceVg and concentration solution space WC are also described as below:
(62)Vg≔uϵH1Ω2∣u=gDon∂ΩD,(63)WC≔ωϵH1Ω2∣ω=CDon∂ΩD.


### 2.7. Weak Formulations

The weak formulations of the model equations are derived for this purpose from Theorem5.Theorem 5.
(Green’s formula) letΩ be an open bounded regular set of class C2. If v and u are functions of H1Ω2, then they satisfy ∫ΩΔuvdx=−∫Ω∇u:∇vdx+∫Ωv∇u.nds, where n is the outward unit normal to ∂Ω and [39]
(64)∇u:∇v=∂u1∂x∂u1∂x+∂u1∂y∂u1∂y+∂u2∂x∂u2∂x+∂u2∂y∂u2∂y=tr∇u∇vT.The weak form of the system equation (52) is as follows.To findu,p,CϵVg×Ψ×WC such that equation (65) is expressed as
(65)∫Ω−αΔu+u.∇uvdx+∫Ωv.∇pdx=∫Ωf.vdx,forallvϵV0,∫Ω∇.uψdx=0,forallψϵΨ,∫Ω∂C∂tωdx−∫ΩDΔCωdx+∫Ωu.∇Cωdx=∫ΩFωdx,forallωϵW0.By applying Green’s formula on the viscous and pressure terms of the momentum balance equation, and on the diffusive term of the convection-diffusion equation, the following is obtained:(66)∫Ωα∇u:∇v+u.∇u.v−p∇.vdx−∫∂ΩNgN.vds=∫Ωf.vdx,forallvϵV0,∫Ω∂C∂tω+D∇C.∇ω+u.∇Cωdx−∫∂ΩNDωCNds=∫ΩFωdx,forallωϵW0.In defining the bilinear forms as in expressions(67)αu,v:V0×V0⟶ℝ,bv,p:V0×Ψ⟶ℝ,bu,ψ:V0×Ψ⟶ℝ,cC,ω:W0×W0⟶ℝ,eu,C,ω:V0×W0×W0⟶ℝ,by
(68)αu,v=∫Ωα∇u:∇vdx,bu,ψ=∫Ωψ∇.udx,bv,p=−∫Ωp∇.vdx,cC,ω=∫Ω∇C:∇ωdx,eu,C,ω=∫Ωu.∇C.ωdx,and the trilinear form
(69)du,u,v:V0×V0×V0⟶ℝ,by
(70)du,u,v=∫Ωu.∇u.vdx.Equation (66) is expressed as
(71)αu,v+du,u,v+bv,p=f,v+gN,v,forallvϵV0,(72)bu,ψ=0,forallψ∈Ψ,(73)C,ω+cC,ω+eu,C,ω=F,ω+Dω,CN,forallωϵW0.The existence and uniqueness of solutions to problems in equations (71) and (73) are briefly discussed.


### 2.8. Existence and Uniqueness of Solutions with the Weak Formulation

The existence and uniqueness of solutions to the weak formulation (equation (71)) is guaranteed by the Lax-Milgram theorem [39]. Alternatively, from the weak formulation (equation (71)), the divergence of free subspace H01Ω2 when introduced can be given by
(74)H≔vϵV0∣∇v=0,and the weak formulation (equation (71)) with the assumption that gN,v=0 can be formulated as follows:
(75)αu,v+du,u,v=f,v,forallv∈H.Lemma 6.
Letu be a solution to a problem in equation (75), the there exists a unique p∈Ψ such that u,p is a solution of problem (equation (71)) [40]. The existence of the solution u in Eequation (75) can be established from the continuity and coercivity arguments which guarantee well-posedness of the problem [40]. For this purpose, the following lemmas and theorems are stated.Lemma 7.
The trilinear formd is continuous.Proof.
The proof follows from the Cauchy-Schwarz inequality, Sobolev embedding theorem, and Holder’s inequality, from which the following can be written asdu,v,w≤cuHvHWH [41].Lemma 8.
The trilinear formd has the following properties:
(a)
du,v,v=0,foru∈V0,v∈H,(b)
du,v,w=−du,v,w,forallu∈V0,v,w∈H. [41]Theorem 9.
There exists a weak solutionu to the Navier-Stokes problem (equation (36)) and a constant c>0 such that
(76)uH≤cαfL2.Proof.
The proof which follows from the Galerkin method involves choosing a finite dimension subspaceV0κ of V0, κ∈N and writing the solution uκ as a linear combination of the functions that form a basis for V0κ since the sequence uκ is bounded in V0 (result from the Galerkin equations), applying the results of Lemma 8, Poincare’s inequality, and the Bolzano Weierstrass theorem, there exists a weakly convergent subsequence uκj such that
(77)uκj⟶uinV0ask⟶∞.
Using compactness and embedding argument, it can be concluded thatu∈V0 is a weak solution to Navier-Stokes problem (equation (36)) [42].The uniqueness of the solutionu is discussed in Theorem 10 below.Theorem 10.
There exists a constantc=cΩ>0, such that the solution of the Navier-Stokes problem (equation (49)) is unique, if α2≥cfL2. Proof: [42]. Remarks on Theorem 10 imply that the uniqueness of the solution to the Navier-Stokes Problem (equation (36)) is guaranteed for sufficiently small data.To establish the existence and uniqueness of solutions to the weak formulation (Equation (73)), Galerkin’s method which involves constructing solutions of certain finite dimensional approximations to equation (73) and then passing to limits can be used [43] in which the procedure involves constructing approximate solutions, deriving energy estimates for the approximate solutions, and then convergence of approximate solutions to a solution.


### 2.9. Numerical Formulations

A numerical model for the groundwater pollute transport problem is derived by discretizing the governing PDE’s by the finite element method (FEM) and in time using the implicit Euler finite difference method.


### 2.10. Approximation with Mixed Finite Element Method

The ability to use a mesh of finite elements to accurately discretize a domain of any size or shape and the possibility of constructing higher order approximating polynomials for greater accuracy makes the finite element method a powerful tool in computational fluid dynamics. FEM is a numerical technique for finding approximate the solutionsux,y,px,y and Cx,y by
(78)ux,y=u0Φ0x,y+u1Φ1x,y+u2Φ2x,y+⋯+ukΦkx,y,px,y=p0ϕ0x,y+p1ϕ1x,y+p2ϕ2x,y+⋯+pkϕkx,y,Cx,y=C0ϕ0x,y+C1ϕ1x,y+C2ϕ2x,y+⋯+Ckϕkx,y,whereΦk and ϕk are shape functions. It also uses the idea that integrals in the weak form (equation (65)) can be broken up into a sum of integrals over an arbitrary collection of disjoint subdomains whose union is the original domain. This implies that the problem can be treated locally and then assembled to the global system.Below is an algorithm for the implementation of the FEM used in the present work.Algorithm 1: For implementing FEM
(a) Obtain the weak formation(b) Discretize the domain in space to obtain the discrete weak formulation(c) Select shape functions (Galerkin’s method)(d) Compute element stiffness, mass, convective, div-grad, and load matrices using a local to global transformation system for shape functions(e) Assemble the global system by adding all the contributions of the local system(f) Implement boundary conditions(g) Solve a global system of equations for the unknownsThe discrete weak formulation is now obtained below.AssumeΩ⊂ℝ2 is polygonal so that it can be tessellated with a set of triangles. We now define a triangulation.Definition 11.
LetΩ⊂ℝ2 be polygonal. A triangular mesh or a triangulation of Ω¯ is a set ℘h of (nondegenerate) N simplices Ki1≤i≤N such that
(79)Ω¯=Ui=1NKi,Ki∈℘h,whereh=maxK∈℘hhk,
andhk is the maximum diameter of triangle K [39] since the spaces H1Ω2,L2Ω, and H1Ω are infinite dimensional; the approximation following finite dimensional subspace is in
(80)Vgh⊂H1Ω2,(81)Ψh⊂L2Ω,(82)Wch⊂H1Ω,which depend on the discretization parameter h. The discrete weak formulation for equation (83) is formulated as follows:
Finduh,phCh∈Vgh×Ψh×Wch,
such that(83)∫Ωα∇uh:∇vh+uh∇uhvh−phdx−∫∂ΩNgNvhdx=∫Ωf.vhdxforallvh∈V0h,(84)∫Ω∂Ch∂tωh+DVCh∇ωh+uh∇Chωhdx−∫∂ΩNDωhCNds=∫ΩFωhdxforallωh∈W0h,(85)∫Ω∇uhψhdx=0forallψh∈Ψh,
Defining bilinear forms(86)auh,vh:V0h×V0h⟶ℝ,buh,ph:V0h×Ψ0h⟶ℝ,buh,ψh:V0h×Ψ0h⟶ℝ,cCh,ωh:W0h×W0h⟶ℝ,euh,Ch,ωh:V0h:W0h×W0h⟶ℝ.
By the expressions in (87)
(87)auh,vh=∫Ωα∇uh:∇vhdx,cCh,ωh=∫ΩD∇Ch:∇ωhdx(88)buh,ph=−∫Ωph∇.vhdx,buh,ψh=∫Ωψh∇:∇.uhdx,(89)euh,Ch,ωh=∫uh∇Ch.ωhdxand the trilinear form
(90)duh,uh,vh:V0h×V0h×V0h⟶ℝ,by
(91)duh,uh,vh=∫Ωuh∇uh.vhdx,
equations (84) can then be expressed in
(92)auh,vh+duh,uh,vh+buh,ph=f,vh+gN,vh,forallvh∈V0h,buh,ψh=0,forallψh∈Ψh,Ch,ωh+cCh,ωh+eub,Ch,ωh=F,ωh+gN,vh,forallωh∈W0h.Ifθjj=1nu is a set of vector-valued-basis functions, and qkk=1np and lmm=1nc scaler basis functions that span the finite dimensional space V0h, ψh and W0h, respectively, such that
(93)θjxi=1,ifi=j,0,ifi≠j,qkxi=1,ifi=k,0,ifi≠k,lmxi=1,ifi=m,0,ifi≠m, then, the discrete velocity uhx, pressure phx, and concentration Chx,t components are written as a linear combination of their respective basis functions as follows:
(94)uhx=∑j=1nuujθjx+∑j=nu+1nu+n∂ujθjx,(95)phx=∑k=1nppkqkx,(96)Chx,t=∑m=1nccmtlmx+∑m=nc+1nc+n∂cmtlmx, where nu,np, and nc are the number of unknowns for velocity, pressure, and concentration, respectively, n∂ is the number of nodal points where velocity or concentration is known on the boundary ∂Ω, and uj and pk are scalar coefficients for the velocity and pressure basis functions, respectively, while Cml is the time-dependent coefficient for the concentration basis function. Substituting the expansions in equations (94) and (95) into equations (84), together with vh=θi and ψh=qk, we obtain equation (97) as
(97)α∑j=1nuuj∫Ω∇θi:∇θjdx+∑j=1nuuj∫Ωθiuh.∇θjdx−∑k=1nppk∫Ωqk∇.θidx=∫Ωf.θidx−∫∂ΩNgN.θids−α∑j=nu+1nu+n∂uj∫Ω∇θi:∇θjdx−∑j=nu+1nu+n∂uj∫Ωθi.uh.∇θjdx, and
(98)−∑j=1nuuj∫Ωqk∇.θjdx=∑j=nu+1nu+n∂uj∫Ωqk∇.θjdx.Then, equations (97) and (98) can be rewritten in
(99)αAu+Buu+YTp=s,Yu=t, or in block matrix form as in
(100)αA+BuYTY0up=st,where(101)A=ai,j,ai,j=∫Ω∇θi:∇θjdx,(102)Y=yk,j,yk,j=−∫Ωqk∇θjdx,(103)Bu=bui,j,bui,j=∫Ωθiuh.∇θjdx, for i,j=1,⋯nu and k=1,⋯,np.The right-hand vectors are expressed in(104)s=si,si=∫Ωf.θidx−∫∂ΩNgN.θids−α∑j=nu+1nu+n∂uj∫Ω∇θi:∇θjdx−∑j=nu+1nu+nθuj∫Ωθi.uh∇θjdx,(105)t=tk,tk=∑j=nu+1nu+n∂uj∫Ωqk∇.θjdx.


### 2.11. Discretization of the Diffusion-Convection Equation

The diffusion-convection equation is discretized in space using finite elements and in time using an implicit finite difference scheme. In this scheme, a system of linear equations is solved to calculate the valuesCk+1jj∈ℤ as functions of the previous values Ckjj∈ℤ at each time step. Substituting the expansion in equation (96) into equation (84) together with wh=lr, the following equation (96) is obtained:
(106)∑m=1nc∂C∂t∫Ωlmlrdx+∑m=1nccm∫ΩD∇lm∇lrdx+∑m=1ncCm∫Ωlruh∇lmdx=∫∂ΩNDCNlrds+∫ΩlrFdx−∑m=nc+1nc+nΩ∂Cm∂t∫Ωlmlrdx−∑m=nc+1nc+nθCm∫ΩD∇lm∇lm∇lrdx−∑m=nc+1nc+nθCm∫Ωlruh∇lmdx.Next, discretize in time by substituting for(107)∂Cm∂t=Ck+1−CkΔt, in equation (106), which gives
(108)∑m=1ncCk+1∫Ωlmlrdx+∑m=1ncDΔtCm∫Ω∇lm∇lrdx+∑m=1ncΔtCm∫Ωlruh∇lmdx=∫∂ΩNΔtDCNlrdx−∫ΩΔtlrFdx−∑m=nc+1nc+n∂Ck+1−Ck∫Ωlmlrdx−∑m=nc+1nc+n∂DΔtCm∫Ω∇lm∇lrdx−∑m=nc+1nc+n∂ΔtCm∫Ωlruh∇lmdx+∑m=1ncCk∫Ωlmlrdx.Equation (108) can be rewritten as shown
(109)EC+DΔtNC+ΔtGuC=z, or in a matrix form as
(110)E+DΔtN+ΔtGuC=z,where in(111)N=nr,m,nr,m=∫Ω∇lm∇lrdx,Gu=gr,m,gr,m=∫Ωlruh∇lmdx,E=er,m,er,mer,m=∫Ωlmlrdx,forr,m=1,⋯..,ncThe right-hand vector in (112) is expressed as
(112)z=zr,zr=∫∂ΩNΔtDCNlrds+∫ΩΔtlrFdx−∑m=nc+1nc+n∂Ck+1−Ck∫Ωlmlrdx−∑m=nc+1nc+n∂ΔtCm∫ΩD∇lmdx−∑m=nc+1nc+n∂Cm∫ΩΔtlruh∇lmdx+∑m=1ncCk∫Ωlmlrdx.The Navier-Stokes system of equation (100) is first solved, and the obtained solution u is utilized to solve the system of the diffusion-convection equation (110). A potential source of instability of solutions to the diffusion-convection equation is now discussed. Methods that are applied to solve problems where no convection is present may totally fail when applied to convection-dominated problems. [44]. For instance, when the diffusivity coefficient D of the diffusion-convection equation is very small, the diffusion-convection equation becomes convective dominated and for such a case, an artificial diffusive term is added through a process known as upwinding to obtain stable solutions.


### 2.12. Streamline Upwind Petrov-Galerkin Method

Solutions to the convection-diffusion equation in equation (66) may develop spurious oscillations (if convection dominated) unless the exact solution happens to be globally smooth [45]. An in-depth analysis on the deficiencies of the classical Galerkin approach in obtaining solutions of convection-dominated transport problems is discussed in [46] The finite element method requires the addition of and an extra term that balances the diffusive terms to obtain stable numerical solutions. A popular and efficient remedy is to augment the Galerkin finite element method formulation of the diffusion-convection equation in equation (83) by an extra term that adds artificial diffusion. This extra term is evaluated over the interior of Ω and is a function of the residual in (113)
(113)RC=∂C∂t−DΔC+u∇C−F,of the PDE in
(114)∂c∂t−DΔC+u.∇C=FinΩT≔0,T×Ω.This consistently stabilized method can take the form of(115)∫Ω∂Ω∂twh+D∇Ch∇wh+uh∇Chwhdx−∫∂ΩNDwhCNds+∑e=1nel∫ΩeRCh/dt℘wh⏟stabilizationterm=∫ΩFwhdx,whereΩe is the interior of Ω, nel the number of interior points in Ω, and RCh is the residual of the diffusion-convection equation (114) and defined in
(116)RCh=∂Ch∂t−DΔCh+uh.∇Ch−F.℘whis an operator applied to the test function wh and
(117)τ=D¯u2,is the stabilization parameter and D¯ defined as the artificial diffusion coefficient in equation (117). The choice of
(118)℘wh=∂wh∂t+u.∇wh, in equation (115), describes what is known as the streamline upwind Petrov-Galerkin (SPUG) finite element method.The word streamline implies that the stabilization parameter is selected in such a way that the artificial diffusion term added is in the direction of flow and not perpendicular to it because convective transport occurs along the streamlines, [40]. The idea of adding diffusion along the flow lines helps to avoid what is called crosswind diffusion [46]. The development of a general theory for selecting D¯ optimally is still an area of research. In equation (115), the stabilization term contains a parameter τ whose value significantly influences the quality of the numerical solution. An optimal choice of τ for convection-dominated equations is usually not known because most of the existing procedures are based on somewhat heuristic procedures, [47]. Detailed discussions on these stabilization schemes are provided [24, 44, 46, 47–54].


### 2.13. Choice of the Finite Element SpacesV0h, Ψh, and W0h

For stability, the velocity is approximated using piecewise quadratic functions while the pressure is estimated using piecewise linear functions [41]. Concentration is also approximated using piecewise linear functions. The discrete spaces are defined in
(119)Vgh≔v∈Vg∩CΩ¯:v∣Ki∈P2KforallKi∈℘h,WCh≔w∈WC∩CΩ¯:wKi∈P1KforallKi∈℘h,and in
(120)V0h≔v∈Vg∩CΩ¯2:v∣Ki∈P2KforallKi∈℘handv=0on∂ΩD,Ψh≔ψ∈Ψ∩CΩ¯:ΨKi∈P1KforallKi∈℘h,W0h≔w∈Wg∩CΩ¯:w∣Ki∈P1KforallKi∈℘handw=0on∂ΩD.Approximating polynomials are chosen in such a way that the velocity is computed at the vertices of the triangle and the midpoints of the edges, while pressure and concentration are computed at the triangle vertices only as shown in Figures4 and 5. The polynomials are built such that they are described on one node and vanish with the others. The quantity of nodes related to a triangular detail is identical to the quantity of nearby ranges of freedom. Convergence houses are ruled with the aid of using the steadiness concerns expressed in the ellipticity requirement of the inf-sup circumstances of Brezzi and Babuska that is a generalization of the coercivity circumstance implied in the Lax-Milgram theorem [41]. For the Navier-Stokes equations, one such preference of a strong and convergent family of finite detail spaces used for approximating velocity and pressure, which can be inf-sup stable is the Taylor-Hood detail p2,p2 used in this study work, [40, 41].Figure 4
Approximate space for pressure and concentration.Figure 5
Approximate space for velocity.


### 2.14. Solving the Discrete Algebraic System

The algebraic system in equation (117) is nonlinear and hence requires an iterative method with a linearized problem being solved at each step. Starting with an appropriately chosen initial guess u0,p0 (the solution of the Stokes flow), the next guess u1,p1 is obtained by solving the resulting linear system. The process is repeated until the difference between the solutions obtained at iteration step N+1 and N is greater than a given value of tolerance. A summary of the above steps can be seen in the algorithm (Figure 6).Figure 6
Steps of the set-up algorithm.Similarly, using an initial guessCk for the diffusion-convection equation, we solve for Ck+1 in
(121)∫ΩCk+1.w−Ck.w+D∇Ck+1.∇wΔt+Δtu.∇Ck+1w−FwΔtdx−∫∂ΩNDw∇Ck+1.nΔtds=0.The above algorithms are implemented in FreeFem and some MATLAB codes in addition are utilized.


## 2.1. A PDE Model for Groundwater Flow

To achieve the objectives of groundwater flow and pollute transport problem, a groundwater flow equation based on an incompressible steady 2-dimensional Navier-Stokes equation and a 2-dimensional convection-diffusion equation describing pollute transport is constructed. The Navier-Stokes equations are derived on the basis of the conservation principles of mass (continuity) and linear momentum (Newton’s second law of motion). The first law of thermodynamics (energy) is not considered because of the assumption that heat transfer in groundwater is negligible. Finally, the boundary conditions both Dirichlet and Neumann are stated, together with reasonable assumptions. The density of water is assumed constant, that is, it is not affected by variations in concentration. The pollute transport equation is derived based on a Fickian model which assumes that the waste is transported mainly due to the concentration gradient and that the dispersive flux occurs in a direction from higher concentrations towards lower concentrations [34]. In this work, the PDE model for pollute transport in groundwater was formulated by first deriving an equation for groundwater flow using Navier-Stokes equations.


## 2.2. Derivation of the Navier-Stokes Equation

The flow of water is described using steady incompressible Navier-Stokes equations. These equations are derived on the basis of the conservation principles of mass (continuity) and linear momentum. LetΩ∈ℝ2 describe the domain covered by the flowing water. A fluid particle X (Figure 3) is considered to be moving through x∈ℝ2 at time t.Figure 3
Mapping from the Langrangian coordinate to Eulerian coordinate.The Lagrangian or material coordinate isX. It is the observed trajectory of a specific fluid particle as the fluid flows while x is a Eulerian coordinate because the fluid particle is observed through a fixed region of interest [35].Also, letωo⊂Ω be a subset of t=0. Define the function ϕ:ℝ2⟶ℝ2 such that det ∇xϕX>0 for all X∈Ω the change of the particle’s position is described by equation (4).
(4)ωt=ϕX,t:X∈ω0=ϕω0,t.Defineu=ux,t to be a vector field on Ω at time t, ρx,t the mass density, ρRX,t the mass density of the reference frame f¯x,t the force density, tx,t,n the surface force density also known as traction or contact forces, and n the outward vector to Ωt. By the conservation of mass property, the total mass remains constant as shown in equation (5).
(5)∫ω0ρRXdx=∫ωtρx,tdx=constantforallt>0.Since the total mass remains constant, the rate of change of mass inωt is given by
(6)ddt∫ωtρx,tdx=0.To simplify the term on the left-hand side of equation (6), Reynolds’ transport (Theorem 1) leading to equation (7) was used.Theorem 1.
For an arbitrary single-valued scalar functionρ=ρx,t with continuous derivates and an arbitrary control region ωt with boundary ∂ωt, outward-pointing unit-normal n, and u the local velocity of the fluid across ∂ωt, the following integral relation holds:
(7)ddt∫ωtρx,tdx=∫ωt∂ρx,t∂tdx+∫∂wtρx,tu.nds,
whereds is an element of length on ∂ωt.Source: Reinstra and Hirschberg [36] and Powers [37].In simple terms, the rate of change of a quantity,ρ=ρx,t, in ωt is a result of the rate at which the quantity is being produced in ωt plus the net amount that crosses the boundary ∂ωt Reynolds’ transport theorem is used to translate the integral conservation laws into differential conservation laws. This helps in formulation of the basic conservation laws of fluid dynamics. To simplify calculations, the boundary integral on the right-hand side of equation (7) is transformed to the domain integral using Gauss divergence (Theorem 2).Theorem 2.
Letsωt a 2-dimensional domain with boundary ∂ωtds an element of length (arc length) on ∂ωt and associated unit outward normal vector n with u a vector field on ωt then ∫ωt∇.udx=∫∂wtu.nds, source: [36]. Gauss’ divergence theorem 1.1 shows the relationship between a domain and surface integral. Applying the Gauss divergence theorem, the right-hand side of equation (7) can be written as
(8)∫∂wtρx,tu.nds=∫ωt∇.ρx,tudx.
Substituting equations (8) into (7) yields equation (9) as
(9)ddt∫ωtρx,tdx=∫ωt∂ρx,t∂tdx+∫ωt∇.ρx,tudx,=∫ωt∂ρx,t∂t+∇.ρx,tudx.
Applying equation (6) to equation (9),
(10)⇒∫ωt∂ρx,t∂t+∇.ρx,tudx=0.Since the integral over ωt is arbitrary, equation (10) can be written as
(11)∂ρx,t∂t+∇.ρx,tu=0Equation (11) is known as the continuity equation. Applying the incompressibility condition ρx,t=ρ=constant to equation (11), the continuity equation now becomes ∇.u=0. The conservation of momentum equation was derived from Newton’s second law of motion. From f¯=ma=d/dt∫ωtρx,tdxu=d/dt∫ωtρx,tudx, where f¯ is the force action on the particle, m is the mass, a is the acceleration, and u is the velocity of the particle. The total force f¯ acting on the subdomain ωt is the sum of the field force fx,t and the surface force tx,t,n.By Newton’s second law of motion,d/dt∫ωtρudx=∫ωtρfx,tdx+∫∂wttx,t,n.nds, where ds is an element of length on ∂ωt. Considering the action of the body or field forces fx,t and the normal component of the surface forces tx,t,n=σx,t.n, where σx,t∈ℝ2×2 is the fluid stress tensor, the formula is written as
(12)ddt∫ωtρudx=∫ωtρfx,tdx+∫∂wtσx,t.nds.Considering theith component of the vectors ux,t,f¯x,t, and σi,∗x,t, where σi,∗x,t denotes the ith row of the matrix σx,t, equation (12) is rewritten as
(13)ddt∫ωtρuidx=∫ωtρfidx+∫∂wtσi,∗.nds,fori=1,2.Using the divergence theorem, the boundary integral on the right-hand side of equation (13) is written as
(14)∫∂wtσi,∗.nds=∫∂wt∇.σi,∗dx.Moreover, by Theorem1, the term on the left-hand side of equation (13) can be written as
(15)ddt∫ωtρuidx=∫wt∂ρui∂t+∇.ρuiudx,=∫ωt∂ρui∂t+ui.∇.ρu+∇ui.ρudx.And by the incompressibility condition, the following is obtained:(16)ddt∫ωtρuidx=∫wt∂ρui∂t+u.∇ρuidx, substituting equations (14) and (16) into equation (13),
(17)∫wt∂ρui∂t+u.∇ρuidx=∫ωtρfidx+∫∂wt∇.σi,∗.dx,⟹∫wt∂ρui∂t+u.∇ρui−ρfi−∇.σi,∗dx=0.Since the integral overωt is arbitrary, equation (17) becomes
(18)∂ρui∂t+u.∇ρui−ρfi−∇.σi,∗=0,where∇.σ=∂σ11∂x+∂σ12∂y∂σ21∂x+∂σ22∂y,∇u=∂u1∂x+∂u1∂y∂u2∂x+∂u2∂y.Introducing constitutive relations into equation (18), σ=−pI+τ and p∈ℝ,I∈ℝ2×2; τ∈ℝ2×2=λ∇.uI+2μD, where λ is the volume viscosity, μ is coefficient of shear viscosity, and we set the deformation tensor D in
(19)D=12∇u+∇uT, is the small strain tensor [38]. Substituting for p, D and σ for ∇.σ in equation (18),
(20)∇.σ=∇.−pI+λ∇.uI+μ∇u+∇uT=−∇p+λ+μ∇∇.u+μ∇.∇u,where∇.∇u=∇2u=Δu is the Laplacian of u. Applying the incompressibility condition to equation (20),
(21)⟹λ+μ∇∇.u=0.Therefore, equation (21) is obtained as
(22)∇.σ=−∇p+μ∇u. Substituting equation (22) into equation (18),
(23)⟹∂ρu∂t+u.∇ρu−ρf+∇p−μΔu=0.Assuming the mass densityρ to be constant, equation (23) can be written as
(24)∂u∂t+u.∇u=f−1ρ∇p+vΔu,wherev=μ/ρ is the kinetic viscosity.To make an analysis of the physical problem when the equation parameters of equation (24) change, scaling of the Navier-Stokes equation (24) was carried out. Normalization of these scales resulted in the formulation of dimensionless groups such as Reynolds’, Froude’s, Euler’s, Weber’s, Prandtl’s, and Mach’s numbers which represent the relative importance of various parts of the Navier-Stokes equations and are also key in determining the flow regimes. Scaling also helped to reduce the number of variables. For instance, the combined effect of both viscosity and density can be determined by one dimensionless variable called Reynolds’ number.


## 2.3. Nondimensionalization of Navier-Stokes Equations

Let(25)u~=uU,t~=tT,x~=xL,T=LU,whereU, L, and T are characteristics speed, length, and time, respectively. From expressions in (25), we obtain
(26)∂u∂t=∂Uu~∂Tt~=UT∂u~∂t~=UTut.~In equation (26), the convective termu.∇ubecomes
(27)u.∇u=Uu~.∇Uu~=U2u~.∇u~=U2u~.∂u~∂x~L=U2u~L.∂u~∂x=U2Lu~.∇~u~.The pressure term becomes(28)1ρ∇p=1ρ∂p∂x=1ρ∂p∂x~L=1ρL∂p∂x~=1ρL∂p∂x~=1ρL∇~p.The viscous termvΔu in (29) becomes
(29)vΔu=v∂2u∂x2=v∂2Uu~∂x~L2=vL2∂2Uu~∂x~2=vUL2∂2u~∂x~2=vUL2Δ~u~,(30)∇.u=∇.u~U=U∇.u,~(31)⟹∇.u=U∂∂x.u~=U∂∂x~L.u~=UL∂∂x~.u~,(32)⟹∇.u=UL∇~.u~=0. Substituting equations (26)–(30) into equation (24), we obtain
(33)⟹UTut~+U2Tu~.∇~u~=f−1ρL∇~p+vUL2Δ~u~.Also letp~=1/ρU2P, then
(34)∇~p=∇~ρU2p~=ρU2∇~p.Introducing a nondimensional source term(35)f~=LU2f,and substituting equations (34) and (35) into equation (33) and further substituting T=L/U and multiplying L/U2, we obtain
(36)u~t+u~.∇~u~=f~−∇~p~+1ReΔ~u~.whereRe=LU/v=ρLU/μ is the Reynolds number. It is a dimensionless ratio of inertial to viscous forces used to determine the flow regime. Dropping the tilde notation and assuming steady state, equations (36) and (28) form
(37)−αΔu+u.∇u+∇p=p.And∇.u=0, where α=1/Re, equation (36) represents the dimensionless steady incompressible Navier-Stokes equations.Pollute transport is due to the combined effect of diffusion and advection in a fluid. The PDE describing it as derived from the principle of conservation of mass using Fick’s law [30]. Underground water flows through the pores (space between soil particles) in the soil. Some soils are more compact than others. The ease with which water flows through the soil depends on the soil porosity denoted by ϵ∈0,1 such that
(38)ϵ=TotalvolumeofvoidsinthesoilTotalvolumeofsoil.


## 2.4. Diffusion-Convection Equation

Consider a porous domainΩ∈ℝ2 and ωt⊂Ωt with Cx,t denoting the density or concentration of the pollute in mass per unit area, Qx,t denotes the flux, that is, mass per unit length per unit time crossing the boundary qx,t the rate at which the pollutant is increasing or reducing in ωt due to the source term. The total mass of pollute in ωt is ∫ωtϵCdx and the amount of pollute that crosses the boundary ∂ωt is ∫ωtQ.nds., where ds is an element of length on ∂ωt and n is a unit outer normal vector to ds. The pollutant crosses the boundary ∂ωt in three ways:
(1)
Advection: this is due to the bulk flow of water molecules where pollute particles are carried along the streamlines. It is given byQα=Cv (concentration X velocity)(2)
Molecular diffusion: this is caused by the random motion and collusions of pollutant particles. This occurs even when the fluid mass is at rest(3)
Mechanical dispersion: this is due to the velocity variations caused by the different flow paths that the groundwater takesThe total increase or decrease of the pollutant inωt due to the source terms is given by ∫ωtqdx. By the mass balance principle, the rate of change of the total mass of the pollutant in ωt equals the net rate of pollutant mass that flows through the boundary ∂ωt plus the net rate of increase or decrease of pollutant mass in ωt; this is summarized as
(39)localchange=advection+diffusion+source.Thus,(40)ddt∫ωtϵCdx=−∫ωtQ.nds+∫ωtqdx.By the divergence theorem, the first term on the right-hand side of equation (40) can be written
(41)−∫ωtQ.nds=−∫ωt∇.Qdx,so that equation (41) becomes
(42)ddt∫ωtϵCdx=−∫ωt∇.Qdx+∫ωtqdx.Applying the Reynolds transport theorem, Gauss divergence theorem, and the incompressibility condition, the derivative term on the left-hand side of equation (42) is written as
(43)ddt∫ωtϵCdx=∫ωt∂∂tϵCdx+∫ωt∇.ϵCvdx,which further simplifies to(44)∫ωt∂∂tϵC+ϵv.∇C+∇.Q−qdx=0.Sinceωt is arbitrary, equation (44) implies
(45)ϵ∂∂tϵC+ϵv.∇C+∇.Q−q=0.The contribution to flux from the effect of diffusion (molecular) and solubility (mechanical dispersion) according to Fick’s law:(46)Qd=−Dx,t∇ϵC,whereDx,t is the diffusion coefficient. The equation means that the flux is proportional to the negative gradient of concentration and the pollution will diffuse from a region of high to low concentration. The total flux is the sum of advective, molecular, and mechanical dispersions [34]. Therefore,
(47)Q=Qα+Qd=−Dx,t∇ϵC+Cv,wherevx is the convection velocity also known as advection velocity, substituting equation (47) into equation (45) and assuming the diffusivity Dx,t is constant, ϵ∂/∂tC+ϵv.∇C=−∇.−D∇ϵC+Cv+q, and by the incompressibility condition ∇.v=0, we obtain
(48)ϵ∂∂tC+ϵv.∇C=ϵDΔC−v.∇C+q,which simplifies to(49)∂C∂t=DΔC−u.∇C+F,whereu=ϵ+1ϵ+1/ϵv is the average velocity vector of groundwater obtained from solving equation (37) and F=q/ϵ is the source term. Combining equations (37) and (49), the complete model equations which are a set of highly nonlinear partial differential is as follows;
(50)−αΔu+u∇u+∇p=f,inΩ,∇.u=0,inΩ,∂C∂t−D∇C+u.∇C=F,inΩT≔0,T×Ω.With initial conditionCx,0=C0 and boundary conditions prescribed by
(51)u=gDon∂ΩD−pn+α∇u.n=gNon∂ΩΝ,where∂ΩD and ∂ΩΝ represent the domain boundary prescribed by Dirichlet and Neumann conditions, respectively, and ∂Ω=∂ΩD∪∂ΩΝ.Also, it is expressed as(52)C=CDon∂ΩDand∂C∂n=CNon∂ΩN,whereC0,CD,CN,gD, and gN are known data to be prescribed later. Dirichletian boundary conditions are called essential because they are imposed explicitly on the solution, that is, on the test space, while Neumann boundary conditions are called natural because they are implicitly defined in the variationally weak formulations in which the model is based on the following assumptions:
(a)
It is known a priority regions where groundwater flows(b)
Groundwater flows at a rate which depends on the permeability of the soil(c)
No chemical reactions occur during the interaction of pollutants with waterThe existence of solutions to equation (50) is discussed below.


## 2.5. Existence of Solutions to the Model Equations

The following definitions are used to describe the functional spaces.Definition 3.
The space of real valued functions that are square integral onΩ with respect to the Lebesgue measure dx is denoted by the relationship
(53)L2Ω≔u:Ω⟶ℝ∣∫Ωu2xdx<∞,(54)L2Ω2≔u:Ω⟶ℝ2∣uiϵL2,i=1,2,
The spaceL2Ω2 is equipped with the norm
(55)uL2Ω2=∫Ωux2dx1/2,and scalar products in [39]
(56)u,vL2Ω2=∫Ωux.vxdx.Definition 4.
The Sobolev spaceH1Ω is denoted by the expressions in [39]
(57)H1Ω≔υ:Ω⟶ℝ∣υϵL2Ω,∇υϵL2Ω2,(58)H1Ω2≔v:Ω⟶ℝ2∣uiϵH2,i=1,2.
This space is equipped with the scalar product in(59)u,vH1Ω2=∫Ωuv+∇vdx,that induces the norm in [39]
(60)uL2Ω2=∫Ωu2+∇u2dx1/2.


## 2.6. Functional Spaces

LetV0,ψ,andW0 be test spaces for velocity, pressure, and pollute concentration, respectively, as described below;
(61)V0≔uϵH1Ω2∣u=0on∂ΩD,ψ≔ψϵL2Ωand,W0≔ωϵH1Ω∣ω=0on∂ΩDThe velocity solution spaceVg and concentration solution space WC are also described as below:
(62)Vg≔uϵH1Ω2∣u=gDon∂ΩD,(63)WC≔ωϵH1Ω2∣ω=CDon∂ΩD.


## 2.7. Weak Formulations

The weak formulations of the model equations are derived for this purpose from Theorem5.Theorem 5.
(Green’s formula) letΩ be an open bounded regular set of class C2. If v and u are functions of H1Ω2, then they satisfy ∫ΩΔuvdx=−∫Ω∇u:∇vdx+∫Ωv∇u.nds, where n is the outward unit normal to ∂Ω and [39]
(64)∇u:∇v=∂u1∂x∂u1∂x+∂u1∂y∂u1∂y+∂u2∂x∂u2∂x+∂u2∂y∂u2∂y=tr∇u∇vT.The weak form of the system equation (52) is as follows.To findu,p,CϵVg×Ψ×WC such that equation (65) is expressed as
(65)∫Ω−αΔu+u.∇uvdx+∫Ωv.∇pdx=∫Ωf.vdx,forallvϵV0,∫Ω∇.uψdx=0,forallψϵΨ,∫Ω∂C∂tωdx−∫ΩDΔCωdx+∫Ωu.∇Cωdx=∫ΩFωdx,forallωϵW0.By applying Green’s formula on the viscous and pressure terms of the momentum balance equation, and on the diffusive term of the convection-diffusion equation, the following is obtained:(66)∫Ωα∇u:∇v+u.∇u.v−p∇.vdx−∫∂ΩNgN.vds=∫Ωf.vdx,forallvϵV0,∫Ω∂C∂tω+D∇C.∇ω+u.∇Cωdx−∫∂ΩNDωCNds=∫ΩFωdx,forallωϵW0.In defining the bilinear forms as in expressions(67)αu,v:V0×V0⟶ℝ,bv,p:V0×Ψ⟶ℝ,bu,ψ:V0×Ψ⟶ℝ,cC,ω:W0×W0⟶ℝ,eu,C,ω:V0×W0×W0⟶ℝ,by
(68)αu,v=∫Ωα∇u:∇vdx,bu,ψ=∫Ωψ∇.udx,bv,p=−∫Ωp∇.vdx,cC,ω=∫Ω∇C:∇ωdx,eu,C,ω=∫Ωu.∇C.ωdx,and the trilinear form
(69)du,u,v:V0×V0×V0⟶ℝ,by
(70)du,u,v=∫Ωu.∇u.vdx.Equation (66) is expressed as
(71)αu,v+du,u,v+bv,p=f,v+gN,v,forallvϵV0,(72)bu,ψ=0,forallψ∈Ψ,(73)C,ω+cC,ω+eu,C,ω=F,ω+Dω,CN,forallωϵW0.The existence and uniqueness of solutions to problems in equations (71) and (73) are briefly discussed.


## 2.8. Existence and Uniqueness of Solutions with the Weak Formulation

The existence and uniqueness of solutions to the weak formulation (equation (71)) is guaranteed by the Lax-Milgram theorem [39]. Alternatively, from the weak formulation (equation (71)), the divergence of free subspace H01Ω2 when introduced can be given by
(74)H≔vϵV0∣∇v=0,and the weak formulation (equation (71)) with the assumption that gN,v=0 can be formulated as follows:
(75)αu,v+du,u,v=f,v,forallv∈H.Lemma 6.
Letu be a solution to a problem in equation (75), the there exists a unique p∈Ψ such that u,p is a solution of problem (equation (71)) [40]. The existence of the solution u in Eequation (75) can be established from the continuity and coercivity arguments which guarantee well-posedness of the problem [40]. For this purpose, the following lemmas and theorems are stated.Lemma 7.
The trilinear formd is continuous.Proof.
The proof follows from the Cauchy-Schwarz inequality, Sobolev embedding theorem, and Holder’s inequality, from which the following can be written asdu,v,w≤cuHvHWH [41].Lemma 8.
The trilinear formd has the following properties:
(a)
du,v,v=0,foru∈V0,v∈H,(b)
du,v,w=−du,v,w,forallu∈V0,v,w∈H. [41]Theorem 9.
There exists a weak solutionu to the Navier-Stokes problem (equation (36)) and a constant c>0 such that
(76)uH≤cαfL2.Proof.
The proof which follows from the Galerkin method involves choosing a finite dimension subspaceV0κ of V0, κ∈N and writing the solution uκ as a linear combination of the functions that form a basis for V0κ since the sequence uκ is bounded in V0 (result from the Galerkin equations), applying the results of Lemma 8, Poincare’s inequality, and the Bolzano Weierstrass theorem, there exists a weakly convergent subsequence uκj such that
(77)uκj⟶uinV0ask⟶∞.
Using compactness and embedding argument, it can be concluded thatu∈V0 is a weak solution to Navier-Stokes problem (equation (36)) [42].The uniqueness of the solutionu is discussed in Theorem 10 below.Theorem 10.
There exists a constantc=cΩ>0, such that the solution of the Navier-Stokes problem (equation (49)) is unique, if α2≥cfL2. Proof: [42]. Remarks on Theorem 10 imply that the uniqueness of the solution to the Navier-Stokes Problem (equation (36)) is guaranteed for sufficiently small data.To establish the existence and uniqueness of solutions to the weak formulation (Equation (73)), Galerkin’s method which involves constructing solutions of certain finite dimensional approximations to equation (73) and then passing to limits can be used [43] in which the procedure involves constructing approximate solutions, deriving energy estimates for the approximate solutions, and then convergence of approximate solutions to a solution.


## 2.9. Numerical Formulations

A numerical model for the groundwater pollute transport problem is derived by discretizing the governing PDE’s by the finite element method (FEM) and in time using the implicit Euler finite difference method.


## 2.10. Approximation with Mixed Finite Element Method

The ability to use a mesh of finite elements to accurately discretize a domain of any size or shape and the possibility of constructing higher order approximating polynomials for greater accuracy makes the finite element method a powerful tool in computational fluid dynamics. FEM is a numerical technique for finding approximate the solutionsux,y,px,y and Cx,y by
(78)ux,y=u0Φ0x,y+u1Φ1x,y+u2Φ2x,y+⋯+ukΦkx,y,px,y=p0ϕ0x,y+p1ϕ1x,y+p2ϕ2x,y+⋯+pkϕkx,y,Cx,y=C0ϕ0x,y+C1ϕ1x,y+C2ϕ2x,y+⋯+Ckϕkx,y,whereΦk and ϕk are shape functions. It also uses the idea that integrals in the weak form (equation (65)) can be broken up into a sum of integrals over an arbitrary collection of disjoint subdomains whose union is the original domain. This implies that the problem can be treated locally and then assembled to the global system.Below is an algorithm for the implementation of the FEM used in the present work.Algorithm 1: For implementing FEM
(a) Obtain the weak formation(b) Discretize the domain in space to obtain the discrete weak formulation(c) Select shape functions (Galerkin’s method)(d) Compute element stiffness, mass, convective, div-grad, and load matrices using a local to global transformation system for shape functions(e) Assemble the global system by adding all the contributions of the local system(f) Implement boundary conditions(g) Solve a global system of equations for the unknownsThe discrete weak formulation is now obtained below.AssumeΩ⊂ℝ2 is polygonal so that it can be tessellated with a set of triangles. We now define a triangulation.Definition 11.
LetΩ⊂ℝ2 be polygonal. A triangular mesh or a triangulation of Ω¯ is a set ℘h of (nondegenerate) N simplices Ki1≤i≤N such that
(79)Ω¯=Ui=1NKi,Ki∈℘h,whereh=maxK∈℘hhk,
andhk is the maximum diameter of triangle K [39] since the spaces H1Ω2,L2Ω, and H1Ω are infinite dimensional; the approximation following finite dimensional subspace is in
(80)Vgh⊂H1Ω2,(81)Ψh⊂L2Ω,(82)Wch⊂H1Ω,which depend on the discretization parameter h. The discrete weak formulation for equation (83) is formulated as follows:
Finduh,phCh∈Vgh×Ψh×Wch,
such that(83)∫Ωα∇uh:∇vh+uh∇uhvh−phdx−∫∂ΩNgNvhdx=∫Ωf.vhdxforallvh∈V0h,(84)∫Ω∂Ch∂tωh+DVCh∇ωh+uh∇Chωhdx−∫∂ΩNDωhCNds=∫ΩFωhdxforallωh∈W0h,(85)∫Ω∇uhψhdx=0forallψh∈Ψh,
Defining bilinear forms(86)auh,vh:V0h×V0h⟶ℝ,buh,ph:V0h×Ψ0h⟶ℝ,buh,ψh:V0h×Ψ0h⟶ℝ,cCh,ωh:W0h×W0h⟶ℝ,euh,Ch,ωh:V0h:W0h×W0h⟶ℝ.
By the expressions in (87)
(87)auh,vh=∫Ωα∇uh:∇vhdx,cCh,ωh=∫ΩD∇Ch:∇ωhdx(88)buh,ph=−∫Ωph∇.vhdx,buh,ψh=∫Ωψh∇:∇.uhdx,(89)euh,Ch,ωh=∫uh∇Ch.ωhdxand the trilinear form
(90)duh,uh,vh:V0h×V0h×V0h⟶ℝ,by
(91)duh,uh,vh=∫Ωuh∇uh.vhdx,
equations (84) can then be expressed in
(92)auh,vh+duh,uh,vh+buh,ph=f,vh+gN,vh,forallvh∈V0h,buh,ψh=0,forallψh∈Ψh,Ch,ωh+cCh,ωh+eub,Ch,ωh=F,ωh+gN,vh,forallωh∈W0h.Ifθjj=1nu is a set of vector-valued-basis functions, and qkk=1np and lmm=1nc scaler basis functions that span the finite dimensional space V0h, ψh and W0h, respectively, such that
(93)θjxi=1,ifi=j,0,ifi≠j,qkxi=1,ifi=k,0,ifi≠k,lmxi=1,ifi=m,0,ifi≠m, then, the discrete velocity uhx, pressure phx, and concentration Chx,t components are written as a linear combination of their respective basis functions as follows:
(94)uhx=∑j=1nuujθjx+∑j=nu+1nu+n∂ujθjx,(95)phx=∑k=1nppkqkx,(96)Chx,t=∑m=1nccmtlmx+∑m=nc+1nc+n∂cmtlmx, where nu,np, and nc are the number of unknowns for velocity, pressure, and concentration, respectively, n∂ is the number of nodal points where velocity or concentration is known on the boundary ∂Ω, and uj and pk are scalar coefficients for the velocity and pressure basis functions, respectively, while Cml is the time-dependent coefficient for the concentration basis function. Substituting the expansions in equations (94) and (95) into equations (84), together with vh=θi and ψh=qk, we obtain equation (97) as
(97)α∑j=1nuuj∫Ω∇θi:∇θjdx+∑j=1nuuj∫Ωθiuh.∇θjdx−∑k=1nppk∫Ωqk∇.θidx=∫Ωf.θidx−∫∂ΩNgN.θids−α∑j=nu+1nu+n∂uj∫Ω∇θi:∇θjdx−∑j=nu+1nu+n∂uj∫Ωθi.uh.∇θjdx, and
(98)−∑j=1nuuj∫Ωqk∇.θjdx=∑j=nu+1nu+n∂uj∫Ωqk∇.θjdx.Then, equations (97) and (98) can be rewritten in
(99)αAu+Buu+YTp=s,Yu=t, or in block matrix form as in
(100)αA+BuYTY0up=st,where(101)A=ai,j,ai,j=∫Ω∇θi:∇θjdx,(102)Y=yk,j,yk,j=−∫Ωqk∇θjdx,(103)Bu=bui,j,bui,j=∫Ωθiuh.∇θjdx, for i,j=1,⋯nu and k=1,⋯,np.The right-hand vectors are expressed in(104)s=si,si=∫Ωf.θidx−∫∂ΩNgN.θids−α∑j=nu+1nu+n∂uj∫Ω∇θi:∇θjdx−∑j=nu+1nu+nθuj∫Ωθi.uh∇θjdx,(105)t=tk,tk=∑j=nu+1nu+n∂uj∫Ωqk∇.θjdx.


## 2.11. Discretization of the Diffusion-Convection Equation

The diffusion-convection equation is discretized in space using finite elements and in time using an implicit finite difference scheme. In this scheme, a system of linear equations is solved to calculate the valuesCk+1jj∈ℤ as functions of the previous values Ckjj∈ℤ at each time step. Substituting the expansion in equation (96) into equation (84) together with wh=lr, the following equation (96) is obtained:
(106)∑m=1nc∂C∂t∫Ωlmlrdx+∑m=1nccm∫ΩD∇lm∇lrdx+∑m=1ncCm∫Ωlruh∇lmdx=∫∂ΩNDCNlrds+∫ΩlrFdx−∑m=nc+1nc+nΩ∂Cm∂t∫Ωlmlrdx−∑m=nc+1nc+nθCm∫ΩD∇lm∇lm∇lrdx−∑m=nc+1nc+nθCm∫Ωlruh∇lmdx.Next, discretize in time by substituting for(107)∂Cm∂t=Ck+1−CkΔt, in equation (106), which gives
(108)∑m=1ncCk+1∫Ωlmlrdx+∑m=1ncDΔtCm∫Ω∇lm∇lrdx+∑m=1ncΔtCm∫Ωlruh∇lmdx=∫∂ΩNΔtDCNlrdx−∫ΩΔtlrFdx−∑m=nc+1nc+n∂Ck+1−Ck∫Ωlmlrdx−∑m=nc+1nc+n∂DΔtCm∫Ω∇lm∇lrdx−∑m=nc+1nc+n∂ΔtCm∫Ωlruh∇lmdx+∑m=1ncCk∫Ωlmlrdx.Equation (108) can be rewritten as shown
(109)EC+DΔtNC+ΔtGuC=z, or in a matrix form as
(110)E+DΔtN+ΔtGuC=z,where in(111)N=nr,m,nr,m=∫Ω∇lm∇lrdx,Gu=gr,m,gr,m=∫Ωlruh∇lmdx,E=er,m,er,mer,m=∫Ωlmlrdx,forr,m=1,⋯..,ncThe right-hand vector in (112) is expressed as
(112)z=zr,zr=∫∂ΩNΔtDCNlrds+∫ΩΔtlrFdx−∑m=nc+1nc+n∂Ck+1−Ck∫Ωlmlrdx−∑m=nc+1nc+n∂ΔtCm∫ΩD∇lmdx−∑m=nc+1nc+n∂Cm∫ΩΔtlruh∇lmdx+∑m=1ncCk∫Ωlmlrdx.The Navier-Stokes system of equation (100) is first solved, and the obtained solution u is utilized to solve the system of the diffusion-convection equation (110). A potential source of instability of solutions to the diffusion-convection equation is now discussed. Methods that are applied to solve problems where no convection is present may totally fail when applied to convection-dominated problems. [44]. For instance, when the diffusivity coefficient D of the diffusion-convection equation is very small, the diffusion-convection equation becomes convective dominated and for such a case, an artificial diffusive term is added through a process known as upwinding to obtain stable solutions.


## 2.12. Streamline Upwind Petrov-Galerkin Method

Solutions to the convection-diffusion equation in equation (66) may develop spurious oscillations (if convection dominated) unless the exact solution happens to be globally smooth [45]. An in-depth analysis on the deficiencies of the classical Galerkin approach in obtaining solutions of convection-dominated transport problems is discussed in [46] The finite element method requires the addition of and an extra term that balances the diffusive terms to obtain stable numerical solutions. A popular and efficient remedy is to augment the Galerkin finite element method formulation of the diffusion-convection equation in equation (83) by an extra term that adds artificial diffusion. This extra term is evaluated over the interior of Ω and is a function of the residual in (113)
(113)RC=∂C∂t−DΔC+u∇C−F,of the PDE in
(114)∂c∂t−DΔC+u.∇C=FinΩT≔0,T×Ω.This consistently stabilized method can take the form of(115)∫Ω∂Ω∂twh+D∇Ch∇wh+uh∇Chwhdx−∫∂ΩNDwhCNds+∑e=1nel∫ΩeRCh/dt℘wh⏟stabilizationterm=∫ΩFwhdx,whereΩe is the interior of Ω, nel the number of interior points in Ω, and RCh is the residual of the diffusion-convection equation (114) and defined in
(116)RCh=∂Ch∂t−DΔCh+uh.∇Ch−F.℘whis an operator applied to the test function wh and
(117)τ=D¯u2,is the stabilization parameter and D¯ defined as the artificial diffusion coefficient in equation (117). The choice of
(118)℘wh=∂wh∂t+u.∇wh, in equation (115), describes what is known as the streamline upwind Petrov-Galerkin (SPUG) finite element method.The word streamline implies that the stabilization parameter is selected in such a way that the artificial diffusion term added is in the direction of flow and not perpendicular to it because convective transport occurs along the streamlines, [40]. The idea of adding diffusion along the flow lines helps to avoid what is called crosswind diffusion [46]. The development of a general theory for selecting D¯ optimally is still an area of research. In equation (115), the stabilization term contains a parameter τ whose value significantly influences the quality of the numerical solution. An optimal choice of τ for convection-dominated equations is usually not known because most of the existing procedures are based on somewhat heuristic procedures, [47]. Detailed discussions on these stabilization schemes are provided [24, 44, 46, 47–54].


## 2.13. Choice of the Finite Element SpacesV0h, Ψh, and W0h

For stability, the velocity is approximated using piecewise quadratic functions while the pressure is estimated using piecewise linear functions [41]. Concentration is also approximated using piecewise linear functions. The discrete spaces are defined in
(119)Vgh≔v∈Vg∩CΩ¯:v∣Ki∈P2KforallKi∈℘h,WCh≔w∈WC∩CΩ¯:wKi∈P1KforallKi∈℘h,and in
(120)V0h≔v∈Vg∩CΩ¯2:v∣Ki∈P2KforallKi∈℘handv=0on∂ΩD,Ψh≔ψ∈Ψ∩CΩ¯:ΨKi∈P1KforallKi∈℘h,W0h≔w∈Wg∩CΩ¯:w∣Ki∈P1KforallKi∈℘handw=0on∂ΩD.Approximating polynomials are chosen in such a way that the velocity is computed at the vertices of the triangle and the midpoints of the edges, while pressure and concentration are computed at the triangle vertices only as shown in Figures4 and 5. The polynomials are built such that they are described on one node and vanish with the others. The quantity of nodes related to a triangular detail is identical to the quantity of nearby ranges of freedom. Convergence houses are ruled with the aid of using the steadiness concerns expressed in the ellipticity requirement of the inf-sup circumstances of Brezzi and Babuska that is a generalization of the coercivity circumstance implied in the Lax-Milgram theorem [41]. For the Navier-Stokes equations, one such preference of a strong and convergent family of finite detail spaces used for approximating velocity and pressure, which can be inf-sup stable is the Taylor-Hood detail p2,p2 used in this study work, [40, 41].Figure 4
Approximate space for pressure and concentration.Figure 5
Approximate space for velocity.


## 2.14. Solving the Discrete Algebraic System

The algebraic system in equation (117) is nonlinear and hence requires an iterative method with a linearized problem being solved at each step. Starting with an appropriately chosen initial guess u0,p0 (the solution of the Stokes flow), the next guess u1,p1 is obtained by solving the resulting linear system. The process is repeated until the difference between the solutions obtained at iteration step N+1 and N is greater than a given value of tolerance. A summary of the above steps can be seen in the algorithm (Figure 6).Figure 6
Steps of the set-up algorithm.Similarly, using an initial guessCk for the diffusion-convection equation, we solve for Ck+1 in
(121)∫ΩCk+1.w−Ck.w+D∇Ck+1.∇wΔt+Δtu.∇Ck+1w−FwΔtdx−∫∂ΩNDw∇Ck+1.nΔtds=0.The above algorithms are implemented in FreeFem and some MATLAB codes in addition are utilized.


## 3. Results and Discussions

Consider a hypothetical factory located next to a water body as shown in Figure7. The factory discharges the pollution on part of the land marked Γin. Groundwater flows in the direction Γup to Γdown, where Γup is the boundary between the water body and the groundwater domain while Γdown is the downstream boundary.Figure 7
Domain for the pollute transport problem.A 2-D cut ofΩ in the x−y plane where x represents the horizontal section of Ω and y represents the vertical section of Ω is presented. Figure 8 shows a discretized domain of size x=0 to x=10 and y=0 to y=5. Circular obstacles are inserted to simulate a real groundwater flow process. The domain is now described in detail as follows:
(a)
Γinis the inlet boundary(b)
Γupis the upstream boundary with a baseline concentration of pollute Cp. Γup is the interface between the water body and the computational domain Ω(c)
Γdownis the downstream boundary that is far away such that it can be assumed that the flux is zero(d)
The flux onΓland is proportional to the difference between the concentration of Γland, Cland and the concentration in Ω, C(e)
Diffusivity coefficientD is a constantFigure 8
Computational domain for the pollute transport problem.The following boundary and initial conditions for the diffusion-convection equation (121) are prescribed in
(122)C=CmonΓin,C=CuponΓup,C=0onΓbottom,D∇C.n=0onΓdown,D∇Cn=βCland−ConΓland,β=0.1,and the following boundary conditions. For velocity u=u1u2 in the Navier-Stokes equations are prescribed in expressions
(123)u1=5,yy,u2=0onΓup,u1,u2=0onΓbottomandΓland∪Γin.In addition, no-slip boundary conditions are prescribed on the obstacles inΩ and no boundary condition in prescribed on Γdown, which implies a homogenous Neumann condition on Γdown. In the next section, the Navier-Stokes and diffusion-convection equations are solved to determine the impact of various flow and concentration parameters.


### 3.1. Velocity and Pressure Profile

For Reynolds’ number,Re=10, the velocity profile shown in Figure 9 shows that the maximum velocity of 15.3 occurs along the line y=2.5, which is zero for y=0 and y=5. This is because of the prescribed inflow velocity, which is parabolic and the no-slip conditions imposed on the boundaries. Data in Figure 9 also shows that the velocity increases in regions where the obstacles are widely spaced. The color bar shows the variations in the magnitude of velocity with the color variations from top to bottom in the order of decreasing magnitude. The pressure distribution shown in Figure 10 shows that pressure decreases in the direction of groundwater flow. It has a maximum of 282.36 at x=0. It was reduced from 253.34 at x=2 to 166.29 at x=5, further reduced to 50.123 at 9.6, and finally to 0 when x=10. This is in agreement with the physical interpretation that a fluid flows from a region of higher pressure to a region of lower pressure. Comparing Figures 9 and 10, it is seen that the velocity is highest in regions of higher pressure and lowest in regions of lower pressure.Figure 9
Velocity profile atRe=10.Figure 10
Pressure distribution atRe=10.


### 3.2. Concentration Profile

Prescribing the following boundary and initial conditions in expressions in(124)Cup=10,Cin=100,Cland=1.For the diffusion-convection equation, the transport equation is solved using the velocity obtained from solving the Navier-Stokes equation forRe=10. The concentration profile for diffusivity coefficient D=0.5 is studied. After 5 time steps, a steady-state solution is reached. Figure 11 shows its concentration profile. The concentration on Γin reduces from 100 to 6.3239 at y=4.3 between x=2.5 and x=6.0. It further reduces to 0.81352 at y=4.1 between x=2.3 and x=6.0. The baseline concentration reduces from 10 on Γup to 6.03239 at x,y=0.62,2.6 and x,y=0.6,1.4. It further reduces to 0.81325 at x,y=1.9,2.6. The concentration elsewhere is 0.Figure 11
Concentration,D=0.5, T=5.


### 3.3. Impact of Increasing Inflow and Baseline Concentrations

Keeping the baseline concentration constant and increasing the inflow concentration, the transport equation is solved for the boundary and initial conditions prescribed by the expressions in(125)Cup=10,Cin=500,Cland=1,andD=0.5. At steady state, the concentration profile is shown in Figure 12.Figure 12
Concentration,D=0.5, Cup=10, Cin=500, T=5.Diffusion of molecules depends entirely on movement from regions of high concentration to regions of low concentration. That is, diffusion occurs down the concentration gradient of that molecule. The greater the difference in concentration, the faster the molecule descends along the concentration gradient. If the difference in concentration is not large, the molecules do not move quickly and the diffusion rate decreases. Figure12 shows that the concentration on Γin, reduced from 500 to 31.524 at y=4.3 between x=2.3 and x=6.0. It was then reduced to 3.9664 at y=4.1 between x=2.3 and x=6.0. The baseline concentration reduces from 10 on Γup, to 3.9664 at x,y=0.8,2.5 and x,y=0.8,1.4. The concentration elsewhere is 0. Figures 12 and 13 show that increasing the inflow concentration while keeping the baseline concentration constant increases the amount of the pollutant diffusing from Γin and reduces the amount of pollution due to convective transport from Γup. Keeping the inflow concentration constant and increasing the baseline concentration, the transport equation is solved for the boundary and initial conditions prescribed by expressions in
(126)Cup=50,Cin=100,Cland=1,and D=0.5. At steady state, the concentration profile is shown in Figure 14.Figure 13
Concentration,D=0.5, Cup=10, Cin=100, T=5.Figure 14
Concentration,D=0.5, Cup=50, Cin=100, T=5.Figure14 shows that the concentration on Γin reduces from 100 to 6.03239 at y=4.3 between x=2.3 and x=6.0. It then reduces to 0.81352 at y=4.1 between x=2.3 and x=6.0. The baseline concentration reduces from 50 on Γup to 11.834 at x,y=1.2,2.5. It then reduces to 0.81352 at x,y=2.6,1.4 and x,y=2.8,2.5. Elsewhere, the concentration is 0. Figures 14 and 15 show that increasing the baseline concentration while keeping the inflow concentration constant increases the amount of pollutant due to convective transport from Γup. The amount of pollutant diffusing from Γin is constant. This is because the velocity of groundwater at the boundary Γin is zero.Figure 15
Concentration,D=0.5, Cup=10, Cin=100, T=5.


### 3.4. Impact of Decreasing the Diffusivity Coefficient

Using the boundary and initial conditions prescribed in equation (124) and Re=10 for the Navier-Stokes equations, the transport equation is solved for different diffusivity coefficients.Figure16 shows the steady state concentration for D=0.5 described earlier. Figure 17 shows the steady state concentration for D=0.1. The concentration on Γin reduces from 100 to 6.5495 at y=4.5 between x=2.3 and x=6.0. It then reduced to 1.0524 at y=4.1 between x=2.3 and x=6.0. The baseline concentration reduces from 10 on Γup to 6.5495 at x,y=0.6,1.5 and x,y=0.6,2.5. It was reduced to 1.0524 at x,y=2.5,2.5. The concentration elsewhere is 0.Figure 16
Concentration,D=0.5, Cup=10, Cin=100, T=5.Figure 17
Concentration,D=0.1, Cup=10,Cin=100,T=5.Figure18 shows the steady state concentration profile for D=0.001. The concentration on Γin reduces from 100 to 3.8776 at y=4.3 between x=2.3 and x=6.0, and the baseline concentration reduces from 10 to Γup to 3.8776 at x,y=1.4,1.5 and x,y=1.5,2.5. The concentration elsewhere is 0. Figure 19 shows the concentration profile for D=0.005, the concentration on Γin reduced from 100 to 1.8236 at y=4.3 between x=2.3 and x=6.0, and the baseline concentration reduces from 10 on Γup to 1.8236 at x,y=2.5,2.5. The concentration elsewhere is 0. From Figures 16–19, it is seen that decreasing the diffusivity coefficient D increases the number of pollutants due to convective transport from the boundary Γup and decreases the amount of pollutant that diffuses from the boundary ΓinFigure 18
Concentrations,D=0.001, Cup=10,Cin=100,T=5.Figure 19
Concentrations,D=0.0005, Cup=10,Cin=100,T=5.


### 3.5. Comparison of Various Scenarios for the Management of Polluted Computational Domains for Pollutant Transport Problems and Their Impact on Groundwater Quality

In this study, we present simulation results for four hypothetical case studies. Figure11 shows the results of the concentration profiles and diffusion coefficients (initial and final concentrations in this scenario) of the aquifer for four cases, including the scenario investigated. The velocity and pressure profile [34] concentration profiles for increasing input and base concentrations, decreasing diffusion coefficient and concentration profile for diffusivity coefficient D=0.5. For the three cases, the movement of contaminants can be expressed as an increase in the input concentration, and while maintaining a constant reference concentration, the number of contaminants diffused from Γin increases, and the number of contaminants due to convective movement from Γup decreases. In addition, it can be seen that the decrease in diffusion coefficient D increases the number of pollutants due to convection motion from the boundary Γup and decreases the number of pollutants diffused from the boundary Γin. Concentrations elsewhere are zero, indicating the effect of boundary conditions on groundwater contamination. The same result as above can be achieved by modifying the hypothetical scenario for the contamination site in terms of different concentrations and diffusion coefficients, probing the location of the contaminated site, and keeping 5-step scenario. For the hypothetical scenario, the results show that the groundwater quality can be controlled by the dividing wall in shallow aquifers, but not in deep aquifers. However, lining contaminated areas is an effective way to prevent groundwater contamination in both shallow and deep aquifers [32]. In addition, the use of linings in polluted areas is superior to other methods of collecting and treating wetland wastewater to protect aquifers from contamination. This method can be used taking into account the conditions of the geotechnical characteristics of the soil, the type of lining, the pollution load in these areas, and the type of element.


## 3.1. Velocity and Pressure Profile

For Reynolds’ number,Re=10, the velocity profile shown in Figure 9 shows that the maximum velocity of 15.3 occurs along the line y=2.5, which is zero for y=0 and y=5. This is because of the prescribed inflow velocity, which is parabolic and the no-slip conditions imposed on the boundaries. Data in Figure 9 also shows that the velocity increases in regions where the obstacles are widely spaced. The color bar shows the variations in the magnitude of velocity with the color variations from top to bottom in the order of decreasing magnitude. The pressure distribution shown in Figure 10 shows that pressure decreases in the direction of groundwater flow. It has a maximum of 282.36 at x=0. It was reduced from 253.34 at x=2 to 166.29 at x=5, further reduced to 50.123 at 9.6, and finally to 0 when x=10. This is in agreement with the physical interpretation that a fluid flows from a region of higher pressure to a region of lower pressure. Comparing Figures 9 and 10, it is seen that the velocity is highest in regions of higher pressure and lowest in regions of lower pressure.Figure 9
Velocity profile atRe=10.Figure 10
Pressure distribution atRe=10.


## 3.2. Concentration Profile

Prescribing the following boundary and initial conditions in expressions in(124)Cup=10,Cin=100,Cland=1.For the diffusion-convection equation, the transport equation is solved using the velocity obtained from solving the Navier-Stokes equation forRe=10. The concentration profile for diffusivity coefficient D=0.5 is studied. After 5 time steps, a steady-state solution is reached. Figure 11 shows its concentration profile. The concentration on Γin reduces from 100 to 6.3239 at y=4.3 between x=2.5 and x=6.0. It further reduces to 0.81352 at y=4.1 between x=2.3 and x=6.0. The baseline concentration reduces from 10 on Γup to 6.03239 at x,y=0.62,2.6 and x,y=0.6,1.4. It further reduces to 0.81325 at x,y=1.9,2.6. The concentration elsewhere is 0.Figure 11
Concentration,D=0.5, T=5.


## 3.3. Impact of Increasing Inflow and Baseline Concentrations

Keeping the baseline concentration constant and increasing the inflow concentration, the transport equation is solved for the boundary and initial conditions prescribed by the expressions in(125)Cup=10,Cin=500,Cland=1,andD=0.5. At steady state, the concentration profile is shown in Figure 12.Figure 12
Concentration,D=0.5, Cup=10, Cin=500, T=5.Diffusion of molecules depends entirely on movement from regions of high concentration to regions of low concentration. That is, diffusion occurs down the concentration gradient of that molecule. The greater the difference in concentration, the faster the molecule descends along the concentration gradient. If the difference in concentration is not large, the molecules do not move quickly and the diffusion rate decreases. Figure12 shows that the concentration on Γin, reduced from 500 to 31.524 at y=4.3 between x=2.3 and x=6.0. It was then reduced to 3.9664 at y=4.1 between x=2.3 and x=6.0. The baseline concentration reduces from 10 on Γup, to 3.9664 at x,y=0.8,2.5 and x,y=0.8,1.4. The concentration elsewhere is 0. Figures 12 and 13 show that increasing the inflow concentration while keeping the baseline concentration constant increases the amount of the pollutant diffusing from Γin and reduces the amount of pollution due to convective transport from Γup. Keeping the inflow concentration constant and increasing the baseline concentration, the transport equation is solved for the boundary and initial conditions prescribed by expressions in
(126)Cup=50,Cin=100,Cland=1,and D=0.5. At steady state, the concentration profile is shown in Figure 14.Figure 13
Concentration,D=0.5, Cup=10, Cin=100, T=5.Figure 14
Concentration,D=0.5, Cup=50, Cin=100, T=5.Figure14 shows that the concentration on Γin reduces from 100 to 6.03239 at y=4.3 between x=2.3 and x=6.0. It then reduces to 0.81352 at y=4.1 between x=2.3 and x=6.0. The baseline concentration reduces from 50 on Γup to 11.834 at x,y=1.2,2.5. It then reduces to 0.81352 at x,y=2.6,1.4 and x,y=2.8,2.5. Elsewhere, the concentration is 0. Figures 14 and 15 show that increasing the baseline concentration while keeping the inflow concentration constant increases the amount of pollutant due to convective transport from Γup. The amount of pollutant diffusing from Γin is constant. This is because the velocity of groundwater at the boundary Γin is zero.Figure 15
Concentration,D=0.5, Cup=10, Cin=100, T=5.


## 3.4. Impact of Decreasing the Diffusivity Coefficient

Using the boundary and initial conditions prescribed in equation (124) and Re=10 for the Navier-Stokes equations, the transport equation is solved for different diffusivity coefficients.Figure16 shows the steady state concentration for D=0.5 described earlier. Figure 17 shows the steady state concentration for D=0.1. The concentration on Γin reduces from 100 to 6.5495 at y=4.5 between x=2.3 and x=6.0. It then reduced to 1.0524 at y=4.1 between x=2.3 and x=6.0. The baseline concentration reduces from 10 on Γup to 6.5495 at x,y=0.6,1.5 and x,y=0.6,2.5. It was reduced to 1.0524 at x,y=2.5,2.5. The concentration elsewhere is 0.Figure 16
Concentration,D=0.5, Cup=10, Cin=100, T=5.Figure 17
Concentration,D=0.1, Cup=10,Cin=100,T=5.Figure18 shows the steady state concentration profile for D=0.001. The concentration on Γin reduces from 100 to 3.8776 at y=4.3 between x=2.3 and x=6.0, and the baseline concentration reduces from 10 to Γup to 3.8776 at x,y=1.4,1.5 and x,y=1.5,2.5. The concentration elsewhere is 0. Figure 19 shows the concentration profile for D=0.005, the concentration on Γin reduced from 100 to 1.8236 at y=4.3 between x=2.3 and x=6.0, and the baseline concentration reduces from 10 on Γup to 1.8236 at x,y=2.5,2.5. The concentration elsewhere is 0. From Figures 16–19, it is seen that decreasing the diffusivity coefficient D increases the number of pollutants due to convective transport from the boundary Γup and decreases the amount of pollutant that diffuses from the boundary ΓinFigure 18
Concentrations,D=0.001, Cup=10,Cin=100,T=5.Figure 19
Concentrations,D=0.0005, Cup=10,Cin=100,T=5.


## 3.5. Comparison of Various Scenarios for the Management of Polluted Computational Domains for Pollutant Transport Problems and Their Impact on Groundwater Quality

In this study, we present simulation results for four hypothetical case studies. Figure11 shows the results of the concentration profiles and diffusion coefficients (initial and final concentrations in this scenario) of the aquifer for four cases, including the scenario investigated. The velocity and pressure profile [34] concentration profiles for increasing input and base concentrations, decreasing diffusion coefficient and concentration profile for diffusivity coefficient D=0.5. For the three cases, the movement of contaminants can be expressed as an increase in the input concentration, and while maintaining a constant reference concentration, the number of contaminants diffused from Γin increases, and the number of contaminants due to convective movement from Γup decreases. In addition, it can be seen that the decrease in diffusion coefficient D increases the number of pollutants due to convection motion from the boundary Γup and decreases the number of pollutants diffused from the boundary Γin. Concentrations elsewhere are zero, indicating the effect of boundary conditions on groundwater contamination. The same result as above can be achieved by modifying the hypothetical scenario for the contamination site in terms of different concentrations and diffusion coefficients, probing the location of the contaminated site, and keeping 5-step scenario. For the hypothetical scenario, the results show that the groundwater quality can be controlled by the dividing wall in shallow aquifers, but not in deep aquifers. However, lining contaminated areas is an effective way to prevent groundwater contamination in both shallow and deep aquifers [32]. In addition, the use of linings in polluted areas is superior to other methods of collecting and treating wetland wastewater to protect aquifers from contamination. This method can be used taking into account the conditions of the geotechnical characteristics of the soil, the type of lining, the pollution load in these areas, and the type of element.


## 4. Conclusions and Recommendations

The idea of a depth-dependent initial filtration coefficient is consistent with the following conceptual model: Given a constant flow rate, the larger particles that are more likely to be trapped are more trapped in the porous medium, so most of the larger particles are retained at the top of the conceptual domain as relatively small particles, with the small particles being less likely to be trapped. This disproportionate absorption leads to an increase in the concentration of fine particles in the migratory particle population and the deeper parts of the porous medium have a greater ability to accommodate the particles, but the colloidal population becomes less sticky as the migration distance increases. Finally, particles that are less likely to be captured are more likely to be found in wastewater by transporting longer filters or passing through the medium. As a result, the initial filtration coefficient decreases along the porous medium. Therefore, groundwater is not entirely harmless, like the simulation results in Figures8–19. An increase in the input and base concentrations, as well as a decrease in the diffusivity coefficient D, increases the extent of groundwater contamination. Building a model that takes into account the processes involved in the transport of soil pollutants is complicated. In this model, for example, soil porosity was treated as a constant to facilitate calculation, analysis, and programming, but it is actually a highly variable quantity. This model can be improved by taking into account a 3D transitory equation from Navier-Stokes, adding a variable coefficient of porosity and diffusivity of the soil, increasing the dimensions of the calculation area and the number of obstacles to obtain a better approximation of the geometric structure of the water underground. Industries or others that dispose of toxic substances, such as chemicals, oil, or pharmaceuticals, must treat them before disposal to reduce the level of groundwater contamination. Recommendations of dumping sites should be located far away from water bodies or direction of groundwater flow. To further check the presented simulation results, experimental studies will be conducted.



---

*Source: 7419502-2022-04-12.xml*