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

Most of the real contaminant problems are de ﬁ ned domains that are geometrically complex and can have di ﬀ erent boundary conditions in di ﬀ erent areas. Therefore, it is usually di ﬃ cult to ﬁ nd a solution analytically, so we use the approximate method to generate an approximate function. One answer to this problem is the ﬁ nite element approach (FEM). This study presents a partial di ﬀ erential 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 di ﬀ usion-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 speci ﬁ c and considered water body with the discharge of factory e ﬄ uents on the ground that ends up reducing the quality of groundwater. First, the ﬂ uid ﬂ ow equation is solved to obtain velocity and pressure pro ﬁ les. Steady-state concentration pro ﬁ les were obtained for various values of di ﬀ usion coe ﬃ cient ( D ), baseline, and input concentrations. The results showed that decreasing the di ﬀ usion coe ﬃ cient D increased the number of pollutants for convective transport and decreased the number of pollutants that di ﬀ used 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.


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][11][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 = 2 or 3. These equations solve velocity vector ðu i ðx, tÞÞ1 ≤ 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][21][22], among others. Henry Darcy discovered that water flow through a porous medium was governed by equation (1).
where Q 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 where by: D v /D t = ð∂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).

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 roundoff 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 onedimensional 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.
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 u h is approximated throughout the solution domain as in Equation (3) where ϕ j are suitable bases or shape functions, coefficients c i 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 diffusionconvection 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

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.

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.
The Lagrangian or material coordinate is X. 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).
Define u = uðx, tÞ to be a vector field on Ω at time t, ρðx, tÞ the mass density, ρ R ðX, tÞ the mass density of the reference frame f ðx, tÞ the force density, tðx, 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). ð Since the total mass remains constant, the rate of change of mass in ω t is given by 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: where ds 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 ∂ω t ds an element of length (arc length) on ∂ω t and associated unit outward normal vector n with u a vector field on  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]). 4 Modelling and Simulation in Engineering u:n ds, 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 ð Substituting equations (8) into (7) yields equation (9) as d dt Applying equation (6) to equation (9), Since the integral over ω t is arbitrary, equation (10) can be written as Equation (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 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 f ðx, tÞ and the surface force tðx, t, nÞ.
By Newton's second law of motion, ðd/dtÞ Ð tðx, t, nÞ:n ds, where ds is an element of length on ∂ω t . Considering the action of the body or field forces f ðx, tÞ and the normal component of the surface forces tðx, t, nÞ = σðx, tÞ:n, where σðx, tÞ ∈ ℝ 2×2 is the fluid stress tensor, the formula is written as Considering the i th component of the vectors uðx, tÞ, f ðx, tÞ, and σ i, * ðx, tÞ, where σ i, * ðx, tÞ denotes the i th row of the matrix σðx, tÞ, equation (12) is rewritten as Using the divergence theorem, the boundary integral on the right-hand side of equation (13) is written as Moreover, by Theorem 1, the term on the left-hand side of equation (13) can be written as d dt And by the incompressibility condition, the following is obtained: substituting equations (14) and (16) into equation (13), Since the integral over ω t is arbitrary, equation (17) becomes Introducing constitutive relations into equation (18), where λ is the volume viscosity, μ is coefficient of shear viscosity, and we set the deformation tensor D in is the small strain tensor [38]. Substituting for p, D and σ for ∇:σ in equation (18), where ∇:∇u = ∇ 2 u = Δu is the Laplacian of u. Applying the incompressibility condition to equation (20), Therefore, equation (21) is obtained as Substituting equation (22) into equation (18), Assuming the mass density ρ to be constant, equation (23) can be written as where v = μ/ρ 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.
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 ϵ = Total volume of voids in the soil Total volume of soil : ð38Þ

Diffusion-Convection Equation.
Consider a porous domain Ω ∈ ℝ 2 and ω t ⊂ Ω t with Cðx, tÞ denoting the density or concentration of the pollute in mass per unit area, Qðx, tÞ denotes the flux, that is, mass per unit length per unit time crossing the boundary qðx, 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 Ð ω t Q:n ds:, 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 by Q α = 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 takes The total increase or decrease of the pollutant in ω t due to the source terms is given by Ð ω t qdx: 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 Thus, d dt By the divergence theorem, the first term on the righthand side of equation (40) can be written so that equation (41) becomes 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 which further simplifies to Since ω t is arbitrary, equation (44) implies The contribution to flux from the effect of diffusion (molecular) and solubility (mechanical dispersion) according to Fick's law: where Dðx, 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, where vðxÞ is the convection velocity also known as advection velocity, substituting equation (47) into equation (45) and assuming the diffusivity Dðx, tÞ is constant, ϵð∂/ ∂tÞC + ϵv:∇C = −∇:ð−D∇ðϵCÞ + CvÞ + q, and by the incompressibility condition ∇:v = 0, we obtain 7 Modelling and Simulation in Engineering which simplifies to where u = ððϵ + 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; in Ω, With initial condition Cðx, 0Þ = C 0 and boundary conditions prescribed by where ∂Ω D and ∂Ω Ν represent the domain boundary prescribed by Dirichlet and Neumann conditions, respectively, and ∂Ω = ∂Ω D ∪ ∂Ω Ν .
Also, it is expressed as where C 0 , C D , C N , g D , and g N 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: The space L 2 ðΩÞ 2 is equipped with the norm and scalar products in [39] u, v This space is equipped with the scalar product in that induces the norm in [39] u 2.6. Functional Spaces. Let V 0 , ψ, and W 0 be test spaces for velocity, pressure, and pollute concentration, respectively, as described below; The velocity solution space V g and concentration solution space W C are also described as below: The weak form of the system equation (52) is as follows.
To find ðu, p, CÞ ϵ ðV g × Ψ × W C Þ such that equation (65) is expressed as 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: In defining the bilinear forms as in expressions by and the trilinear form The existence and uniqueness of solutions to problems in equations (71) and (73) are briefly discussed.

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 H 1 0 ðΩÞ 2 when introduced can be given by and the weak formulation (equation (71)) with the assumption that ðg N , vÞ = 0 can be formulated as follows: Lemma 6. Let u 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 form d 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 as dðu, v, wÞ ≤ ck ukHkvkHkWkH [41]. Theorem 9. There exists a weak solution u to the Navier-Stokes problem (equation (36)) and a constant c > 0 such that Proof. The proof which follows from the Galerkin method involves choosing a finite dimension subspace V 0κ of V 0 , κ ∈ N and writing the solution u κ as a linear combination of the functions that form a basis for V 0κ since the sequence ðu κ Þ is bounded in V 0 (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 Using compactness and embedding argument, it can be concluded that u ∈ V 0 is a weak solution to Navier-Stokes problem (equation (36)) [42].
The uniqueness of the solution u is discussed in Theorem 10 below.

Theorem 10.
There exists a constant c = cðΩÞ > 0, such that the solution of the Navier-Stokes problem (equation (49)) is unique, if α 2 ≥ ckf kL 2 . 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.

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.

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 solutions uðx, yÞ, pðx, yÞ and Cðx, yÞ by 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.
The 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.
and h k is the maximum diameter of triangle K [39] since the spaces H 1 ðΩÞ 2 , L 2 ðΩÞ, and H 1 ðΩÞ are infinite dimensional; the approximation following finite dimensional subspace is in which depend on the discretization parameter h. The discrete weak formulation for equation (83) is formulated as follows: 10

Modelling and Simulation in Engineering
Defining bilinear forms By the expressions in (87) and the trilinear form by equations (84) can then be expressed in If fθ j g n u j=1 is a set of vector-valued-basis functions, and fq k g where n u , n p , and n c 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 u j and p k are scalar coefficients for the velocity and pressure basis functions, respectively, while C m ðlÞ is the time-dependent coefficient for the concentration basis function. Substituting the expansions in equations (94) and (95) into equations (84), together with v h = θ i and ψ h = q k , we obtain equation (97) as 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 unknowns and Then, equations (97) and (98) can be rewritten in or in block matrix form as in where for i, j = 1, ⋯n u and k = 1, ⋯, n p . The right-hand vectors are expressed in

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 values ðC Next, discretize in time by substituting for in equation (106), which gives Equation (108) can be rewritten as shown or in a matrix form as where in The right-hand vector in (112) is expressed as 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 diffusionconvection 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.

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 convectiondominated 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 diffusionconvection 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) of the PDE in This consistently stabilized method can take the form of where Ω e is the interior of Ω, n el the number of interior points in Ω, and RðC h Þ is the residual of the diffusionconvection equation (114) and defined in ℘ðw h Þis an operator applied to the test function w h and is the stabilization parameter and D defined as the artificial diffusion coefficient in equation (117). The choice of 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][48][49][50][51][52][53][54].
2.13. Choice of the Finite Element Spaces V h 0 , Ψ h , and W h 0 . 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 and in Approximating polynomials are chosen in such a way that the velocity is computed at the vertices of the triangle 13 Modelling and Simulation in Engineering and the midpoints of the edges, while pressure and concentration are computed at the triangle vertices only as shown in Figures 4 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 ðp 2 , p 2 Þ used in this study work, [40,41].
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 ðu 0 , p 0 Þ (the solution of the Stokes flow), the next guess ðu 1 , p 1 Þ 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).
Similarly, using an initial guess C k for the diffusionconvection equation, we solve for C k+1 in The above algorithms are implemented in FreeFem and some MATLAB codes in addition are utilized.

Results and Discussions
Consider a hypothetical factory located next to a water body as shown in Figure 7. 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.
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: 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.   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 pre-scribed 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.

Concentration Profile. Prescribing the following boundary and initial conditions in expressions in
C up = 10, C in = 100, C land = 1: For the diffusion-convection equation, the transport equation is solved using the velocity obtained from solving the Navier-Stokes equation for Re = 10. The concentration profile for diffusivity coefficient D = 0:5 is studied. After 5 time steps, a steady-state solution is reached. Figure 11 Set u=u 0 and p=p 0 N= 0 Solve for (u N+1 , p=p N+1 ) in

No
Is error less than TOL?  and D = 0:5. At steady state, the concentration profile is shown in Figure 12.
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. Figure 12 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 C up = 50, C in = 100,     and D = 0:5. At steady state, the concentration profile is shown in Figure 14. Figure 14 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.

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. Figure 11 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.

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 . 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.

Data Availability
No raw data were used in support of this study

18
Modelling and Simulation in Engineering