# Numerical Simulation of Pollutant Transport in Fractured Vuggy Porous Karstic Aquifers

Authors: Xiaolin Fan; Shuyu Sun; Wei Wei; Jisheng Kou
Journal: Journal of Applied Mathematics
 (2011)

---

## Abstract

This paper begins with presenting a mathematical model for contaminant transport in the fractured vuggy porous media of a species of contaminant (PCP). Two phases are numerically simulated for a process of contaminant and clean water infiltrated in the fractured vuggy porous media by coupling mixed finite element (MFE) method and finite volume method (FVM), both of which are locally conservative, to approximate the model. A hybrid mixed finite element (HMFE) method is applied to approximate the velocity field for the model. The convection and diffusion terms are approached by FVM and the standard MFE, respectively. The pressure distribution and temporary evolution of the concentration profiles are obtained for two phases. The average effluent concentration on the outflow boundary is obtained at different time and shows some different features from the matrix porous media. The temporal multiscale phenomena of the effluent concentration on the outlet are observed. The results show how the different distribution of the vugs and the fractures impacts on the contaminant transport and the effluent concentration on the outlet. This paper sheds light on certain features of karstic groundwater are obtained.


---

## Body

## 1. Introduction

As is well known, Karst topography is characterized by subterranean limestone caverns, which are carved by groundwater flow. A schematic presentation of a karstic area is showed in Figure1.Figure 1
A schematic presentation of water supply and waste disposal in a karstic area (1) ponor (swallow hole), (2) water well, (3) disposal of municipal waste water into karst aquifer, (4) disposal of industrial waste water into a ponor, (5) solid waste disposal into a karstic depression, (6) karst spring (with limited watershed area), (7) base of karstification (impervious basement), (8) industrial waste water into karst aquifers, (9) municipal waste water into karst aquifers.Karstic area is significantly rich in water resource, the percentage of which for human drinking is up to about 25 and is to be increasing to 50 in the future in the world [1, 2]. Karst formations are cavernous and therefore have very high permeability, resulting in some different features in comparison with non-karstic aquifers. The fluid flow and transport of pollutants occurring in a karstic region have some different features from those in non-karstic areas. With rapid development of local economics and increasing population in the Karst region, numerous pollutants are occurring in some Karst region mainly comprising the following four aspects grouped as municipal, industrial, agricultural, and farraginous. The groundwater in the karst area is seriously affected by human activities on the quality and quantity. Consequently, transport of a contaminant in the Karst topography is an important aspect to be investigated. Although groundwater flow in the Karst terrane has been studied many years in some specific, there are still a number of topics to be researched because of its importance and complexity. One of the most important aspects is the simulation of fluid flow and the pollutant transport in the karstic field.A large number of models have been created to predict water levels and spring discharge, which can reflect some properties of the groundwater flow in the Karst aquifers. Black box model is the simplest one in which no spatial information is included, but can predict spring discharge or other aquifer properties. Time moment analysis is used to relate a time series of inputs (recharge) to a series of outputs [3]. Zaltsberg applied the simple regression models to predict water levels in Karst aquifers [4]. Lack of predictive power is the obvious limitation of these types of models [5, 6].Some deterministic models [5, 6] for groundwater flow are physically derived, and both distributed and lumped parameters are contained therein. The spatial dimension is not contained in lumped models [7] so that they need to solve the only ordinary linear differential equations. Therefore fewer data requirements for parameterization and calibration are provided to simulate a given lumped system than that of distributed counterparts. Quinn et al. [5] and Barrett and Charbeneau [6], for example, developed a system including lumped parameters in groundwater applications generally with single-cell models. Yurtsever and Payne [8] employed a series of linear reservoirs to model Karst aquifers. Normally many researchers apply distributed parameter models to increase the accuracy of predictions or to achieve a high degree of spatial resolution. Also some distributed parameter models using a single equivalent porous medium have been created for the San Antonio portion of Edwards [9] although none performs better than the nine-cell model developed by Wannakule and  Robert Anaya. One of the most sophisticated distributed parameter models was created by Kuniansky and Holligan [10] at the US Geological Survey, which is a finite element model of the Edwards/Trinity aquifer system, in which over 7000 elements were included. Though the high degree of spatial resolution is achieved, difficulties in generating input parameters have limited its use.Some researchers have also developed dual porosity distributed parameter models for Karst system [11]. Conduit and diffuse flow were generally described as separate systems linked by a transfer function in these models. They have the advantage of being able to represent the fast transit and slow depletion often exhibited by Karst aquifers, but at the cost of more than doubling the number of parameters required for calibration.These models demonstrate the evolution in model complexity associated with attempts to increase the accuracy of predictions. The general tendency of the research is to increase the number of cells in thex-y plane while neglecting improvement of which might be achieved by incorporating variation in the vertical direction. This approach has not been consistently successful. The spatially detailed models have been difficult to calibrate and verify. In addition, input data must be developed for each cell; consequently, these models are not used to any great extent by regulatory agencies or other groups. A model developed by Barrett and Charbeneau  [12] is simple to calibrate and use and yet achieves a high degree of accuracy. His modeling effort differs from preceding studies by retaining a simple spatial description of the aquifer, but allowing vertical variations in aquifer properties such as specific yield within cells.Many other numerical models of groundwater flow use finite-difference methods to discretize and solve the governing PDEs [13, 14]. These models often require a highly refined finite-difference mesh to achieve accurate solutions in the areas of interest where gradients vary rapidly in space. Use of a fine mesh over the entire domain can be computationally intensive, and in some cases intractable, while using a variably spaced mesh can lead to cells with a large aspect ratio and refinement in areas where such detail is not needed. In addition, a fine discretization is often needed within models that have already been constructed, and redesigning the entire grid is not feasible. One solution to these predicaments is to use local grid refinement in which the mesh is only refined locally in the area of interest [15, 16].  Mehl and Hill [15, 16] proposed a new method for locally refining block-centered finite-difference grids using iteratively coupled shared nodes (a new method of interpolation) in the context of groundwater flow modeling. His method couples the grids by sharing nodes and iteratively updating the right-hand side of the matrix equations to ensure that heads and fluxes are consistent between both grids. The iterative method presented in their work is a compromise between the accuracy of a variably spaced grid [17, 18] and the speed of the traditional grid refinement methods. In the methods presented in their work, they made use of the Darcy-weighted interpolation to deal with the interface between parent grids and child grids. Two years later [15, 16], they extended the methods to three dimensions mainly by proposing a new 3D interpolation scheme.The selection of the appropriate model to achieve the goals of simulating the groundwater flow in a Karst system is a major task [19]. An appropriate conceptual model should be sufficiently simple so as to be amenable to mathematical treatment, but it should not be too simple so as to exclude those features which are of interest to the investigation at hand. The information should be available for calibrating the model, and the model should be the most economic one for solving the problem at hand.Modeling a system requires a very detailed knowledge of the physical properties and the processes governing water movement. The virtue of a model rests in its ability to predict a general system from incomplete or partial data [20]. The parsimonious model simplifies the representation of the physical structure or of the processes involved [6]. This is especially appropriate in the light of the extraordinary heterogeneity exhibited by Karst aquifers.Recent studies on vuggy fractured porous media simulation are closer to our research.  Fard and Firoozabadi [21] did some comprehensive study of immiscible fluid flow in fractured vuggy porous media. Hoteit and Firoozabadi  [22] showed that multi-component compressible flow in discrete fractured vuggy media can be modeled. Sonnenthal et al. [23] made use of the TOUGHREACT reactive transport software to model the coupled heat transfer and reactive transport processes in porous media. A Crank-Nicolson-Galerkin finite element model is presented in Hossain’s research to simulate nonlinearly the macrophase contaminant transport combined with microphase contaminant transport [24]. Alajmi and Grader [25] gave the relationship between fractured vuggy matrix environment with multiphase flow.  Sun and Geiser [26] and Sun and Wheeler  [27] employed multiscale discontinuous Galerkin and operator-splitting methods to model subsurface flow and transport with anisotropic and dynamic mesh. Sun et al. [28–30] also used the compatible algorithms for coupled flow and transport to simulate the groundwater flow and contaminant transport. Among these approaches the researchers assumed the pressure in a fracture or vuggy element was equal to the pressure in the ambient matrix elements, by using the cross-flow equilibrium.Vugs and fractures are the most important factors in the fluid flow and contaminant transport in the karstic aquifers. We calculate the contaminant transport through the vuggy fractured matrix by using a bulk material property by the distribution coefficientKd, which describes the distribution of contaminant between the liquid and solid phases.  In this paper, the adsorption term is considered to be linear with concentration of the contaminant in the matrix media.In this paper, a karstic aquifers of some size(0,0.80)×(0,0.80) are considered in which the groundwater is clean at first. At some time, a species of contaminant (PCP, i.e., pentachlorophenol) is abruptly infiltrated into the inflow boundary, which lasted T1 hours, and then the clean water is infiltrated into the aquifers again. Among the three phases, the same models of fluid flow and contaminant transport are presented with different boundary conditions and initial conditions which will be described in the sequent sections in this paper. We apply a discrete-fracture model and a discrete-vug model to describe flow and transport processes in fractured porous media and vuggy porous media, respectively. Other than classical discrete fracture or discrete vug model where the cracks or large caves are presented by (n-1) dimensional elements, the discrete fracture or vug model in this paper still uses physically sense n dimensional elements. An adaptive mesh is generated based on this type of model. Then we employ the mixed finite element (MFE) method to approach the second-order partial derivative terms of the flow and transport equations, and an upwind finite volume method (FVM) is used to approximate the convection term in the transport.This paper is organized as follows. First, we mainly give the mathematical models approximately describing the fluid flow and the contaminant diffusion and transport in the fractured vuggy porous media for karst aquifers.  Secondly, the numerical algorithms for the mathematical models are discussed.  We describe in detail the MFE method for the flow equation, the MFE-FVM method for the reactive transport system in the models together with the numerical discretization in time. Then we give some cases in vuggy porous media with vug being differently distributed. For each case, we provide and discuss simulated concentration profiles at different time, as well as velocity and pressure fields. At last, we numerically give the different time scale for concentration in order to analyze the influence on the concentration, pressure, and velocity by vugs and fractures.


## 2. Mathematical Model

Two coupled differential systems are made up of the model equations for the transport of contaminant through the vuggy fractured porous media, that is, the flow equation of the fluid and the transport equations of the contaminant.


### 2.1. Flow Equation

We here consider the flow in vuggy fractured porous media in two-dimension steady flow of single phase, the flow equation of which is obtained from the conservation of total fluid volume and Darcy’s law mathematically given by(2.1)u+K∇p=0,∇⋅u=q,x∈Ω.
Here K denotes the hydraulic conductivity (i.e., the proportionality constant for the flow of water through a porous media, m/s) as follows:(2.2)K=κ⋅gν=κ⋅g⋅ρη.In the two equations above,g is the gravity acceleration, ν and η are the kinematic (s/m2) and dynamic viscosities of the fluid (kg/(m·s)), respectively; ρ and κ are the density of the fluid (kg/m3) and the (absolute) permeability (m2) of the porous media, respectively. The unknowns are p (the pressure head of the fluid mixture (N/m2) and u (the Darcy velocity of the mixture, i.e., the discharge of fluid flowing per unit area, with units of length per time, m/s). We assume that K is uniformly symmetric positive definite and bounded. For simplicity of computation, the computational domain of interest Ω is assumed to be polygonal and bounded in R2 with boundary Γ=∂Ω=ΓD∪ΓN (ΓD, ΓN stand for the Dirichlet boundary and Neumann boundary, resp.). We take the boundary conditions as follows:(2.3)p=pB,x∈ΓD,u⋅n=uB,x∈ΓN,
where pB and uB are the given pressure on ΓD and the normal velocity component on ΓN (pB and uB can be measured in practice in some points on boundary), respectively, and n is the outward norm vector towardsu.


### 2.2. Transport System

The concentration of a contaminant species of interest in the fluid and in the solid together with their relation is described by the following system, which is obtained from the mass conservation of the considered contaminant species:(2.4)∂(ϕcf)∂t+∇⋅(ucf-D(u)∇cf)=r(cf,cs),∂(ϕcs)∂t=-r(cf,cs),x∈Ω,t∈(0,T],cs=Kdcf.
In the system (2.4), ϕ and T stand for the effective porosity and the final simulation time, respectively; ϕ is assumed to be time dependent and uniformly bounded above and below by positive numbers. Kd represents a parameter of the distribution coefficient of the considered contaminant species between the matrix and fluid. D(u), the dispersion-diffusion tensor, contributes from molecular diffusion and mechanical dispersion, and it can be computed by(2.5)D(u)=ϕ{dmI+|u|[dlE(u)+dt(I-E(u))]},
where E(u) denotes the projected tensor onto the u direction, which is calculated by
(2.6)E(u)=1(u12+u22)1/2[u12u1u2u2u1u22],
where |u|=(u12+u22)1/2 represents the Euclidean norm of u=(u1,u2)∈R2; dm (strictly positive) is the molecular diffusivity, dt and dl are the transverse dispersivity and the longitudinal dispersivity, respectively, and both are nonnegative. The unknowns cf as well as cs are the concentrations of the species considered within the fluid and solid, respectively.By summing of the first two equations using the third equation in (2.4), we can obtain(2.7)∂(ϕecf)∂t+∇⋅(ucf-D(u)∇cf)=0,
where ϕe=ϕ+ρKd is calculated for the matrix and for the vugs, respectively, as follows:-in the vugs and/or fractures:ϕ=1, Kd=0, ϕe=ϕ=1.0,-in the matrix:ϕ=ϕm, ϕe=ϕ+ρKd.The boundary conditions for transport system are given by(2.8)(uc-D(u)∇cf)⋅n=cBu⋅n,x∈Γinflow,D(u)∇cf⋅n=0,x∈Γoutflow,
where cB is the prescribed given concentration on the boundary; Γinflow denotes the inflow boundary, and Γoutflow denotes the outflow boundary, defined by(2.9)Γinflow={x∈∂Ω:u⋅n∂Ω≤0},Γoutflow={x∈∂Ω:u⋅n∂Ω>0}.
The initial condition is taken as(2.10)cf(x,0)=c0(x),x∈Ω,
where c0(x) is a prescribed given function for the concentration at t=0, which can be measured for some points in the considered domain in practice.


## 2.1. Flow Equation

We here consider the flow in vuggy fractured porous media in two-dimension steady flow of single phase, the flow equation of which is obtained from the conservation of total fluid volume and Darcy’s law mathematically given by(2.1)u+K∇p=0,∇⋅u=q,x∈Ω.
Here K denotes the hydraulic conductivity (i.e., the proportionality constant for the flow of water through a porous media, m/s) as follows:(2.2)K=κ⋅gν=κ⋅g⋅ρη.In the two equations above,g is the gravity acceleration, ν and η are the kinematic (s/m2) and dynamic viscosities of the fluid (kg/(m·s)), respectively; ρ and κ are the density of the fluid (kg/m3) and the (absolute) permeability (m2) of the porous media, respectively. The unknowns are p (the pressure head of the fluid mixture (N/m2) and u (the Darcy velocity of the mixture, i.e., the discharge of fluid flowing per unit area, with units of length per time, m/s). We assume that K is uniformly symmetric positive definite and bounded. For simplicity of computation, the computational domain of interest Ω is assumed to be polygonal and bounded in R2 with boundary Γ=∂Ω=ΓD∪ΓN (ΓD, ΓN stand for the Dirichlet boundary and Neumann boundary, resp.). We take the boundary conditions as follows:(2.3)p=pB,x∈ΓD,u⋅n=uB,x∈ΓN,
where pB and uB are the given pressure on ΓD and the normal velocity component on ΓN (pB and uB can be measured in practice in some points on boundary), respectively, and n is the outward norm vector towardsu.


## 2.2. Transport System

The concentration of a contaminant species of interest in the fluid and in the solid together with their relation is described by the following system, which is obtained from the mass conservation of the considered contaminant species:(2.4)∂(ϕcf)∂t+∇⋅(ucf-D(u)∇cf)=r(cf,cs),∂(ϕcs)∂t=-r(cf,cs),x∈Ω,t∈(0,T],cs=Kdcf.
In the system (2.4), ϕ and T stand for the effective porosity and the final simulation time, respectively; ϕ is assumed to be time dependent and uniformly bounded above and below by positive numbers. Kd represents a parameter of the distribution coefficient of the considered contaminant species between the matrix and fluid. D(u), the dispersion-diffusion tensor, contributes from molecular diffusion and mechanical dispersion, and it can be computed by(2.5)D(u)=ϕ{dmI+|u|[dlE(u)+dt(I-E(u))]},
where E(u) denotes the projected tensor onto the u direction, which is calculated by
(2.6)E(u)=1(u12+u22)1/2[u12u1u2u2u1u22],
where |u|=(u12+u22)1/2 represents the Euclidean norm of u=(u1,u2)∈R2; dm (strictly positive) is the molecular diffusivity, dt and dl are the transverse dispersivity and the longitudinal dispersivity, respectively, and both are nonnegative. The unknowns cf as well as cs are the concentrations of the species considered within the fluid and solid, respectively.By summing of the first two equations using the third equation in (2.4), we can obtain(2.7)∂(ϕecf)∂t+∇⋅(ucf-D(u)∇cf)=0,
where ϕe=ϕ+ρKd is calculated for the matrix and for the vugs, respectively, as follows:-in the vugs and/or fractures:ϕ=1, Kd=0, ϕe=ϕ=1.0,-in the matrix:ϕ=ϕm, ϕe=ϕ+ρKd.The boundary conditions for transport system are given by(2.8)(uc-D(u)∇cf)⋅n=cBu⋅n,x∈Γinflow,D(u)∇cf⋅n=0,x∈Γoutflow,
where cB is the prescribed given concentration on the boundary; Γinflow denotes the inflow boundary, and Γoutflow denotes the outflow boundary, defined by(2.9)Γinflow={x∈∂Ω:u⋅n∂Ω≤0},Γoutflow={x∈∂Ω:u⋅n∂Ω>0}.
The initial condition is taken as(2.10)cf(x,0)=c0(x),x∈Ω,
where c0(x) is a prescribed given function for the concentration at t=0, which can be measured for some points in the considered domain in practice.


## 3. Numerical Algorithms

In this paper, the system is made up of two parts: the flow equations including the Darcy velocity and the pressure, and the contaminant species transport equations for presenting the evolution of concentration. A triangular mesh for spatial partitioning is employed for accurately approaching the velocity and the concentration around the vugs and fractures. Vugs and fractures are initially characterized by ellipses of different sizes and rectangles of different sizes, respectively. A mixed finite method (FEM) is used to solve the flow equation, based on the triangular mesh; then we solve the transport system by combining the finite volume method (FVM) and MFE method semi-implicitly (implicitly for diffusion and adsorption and explicitly for convection) in time,  getting  ODEs for concentration over time.


### 3.1. Triangular Mesh Generation for Vuggy Fractured Porous Media

Similar to any other numerical modeling, the first step of our algorithm is creating a mesh. The vugs and fractures may be distributed randomly in the domain, so a triangular mesh is required to fit the vuggy fractured porous media much better than a rectangular mesh (though it is easier to generate). A large number of triangular mesh programs are available. Three of them (including DISTMESH, MESHGEN, and TRIANGLE) have much higher mesh quality than any other one in generating adaptive triangular mesh [31]. In this paper, we apply the MESHGEN [32] to generate an adaptive triangular mesh. The fractures are presented in rectangles distributed randomly in the domain considered, and the vugs are described in ellipses distributed randomly.


### 3.2. Mixed Finite Element Method

We employ a mixed finite element (MFE) method for the treatment of the flow equation. MFE method [33] is based on the variational principle that expresses an equilibrium that can be satisfied locally on each finite element. The MFE formulation for the flow equation contains solution for both the scalar variable (pressure) and flux vector (total velocity). Choosing MFE method to approach spaces can satisfy three important properties: local mass conservation, flux continuity, and the same order of convergence (and in some cases super convergence) for both the scalar variable and the flux [33]. MFE method can directly accommodate full permeability tensors and it is more accurate in flux calculation than the conventional finite volume and finite element methods.We first give some symbol notations.  Let(·,·) represent the  L2(Ω) inner product over a domain Ω∈R2 for scalar functions, and (L2(Ω))2 denotes inner product for vector functions. Some necessary spaces are given by(3.1)V=L2(Ω),W=H(div⁡;Ω)={w∈(L2(Ω))2:∇⋅w∈L2(Ω)},W0={w∈H(div⁡;Ω):w⋅n=0on∂Ω},WN0={w∈H(div⁡;Ω):w⋅n=0onΓN}.Based on these spaces, we give the weak formulations for flow equation and reactive equation as follows, respectively.


### 3.3. MFE for Flow Equation

#### 3.3.1. Weak Formulation

For the flow equation, the weak formulation is to findp∈Vand u∈WN0+E(uB) such that (3.2)(K-1u,w)-(∇⋅w,p)=-∫ΓpBw⋅nds,∀w∈WN0,t∈(0,T],(∇⋅u,v)=(q,v),∀v∈V,t∈(0,T],
where E(uB) denotes the velocity extension which lets its component agrees with uB on ΓN. There is the derivation of (3.2) in Appendix A.


#### 3.3.2. MFE Scheme

We apply the RT (Raviart-Thomas) finite element space [34] to approximate the Darcy velocity. As is well known, the rth order RT space for the two-dimensional triangular element is given by(3.3)Vh(K)=Pr(K),Wh(K)=(x1,Pr(K)x2Pr(K))⊕(Pr(K))2,
where ⊕ represents the direct sum.  We employ the case r=0 in this paper. Therefore, Wh(K) has the form(3.4)Wh(K)={w:w=(aK+bKx1,cK+bKx2),aK,bK,cK∈R}.
Restricted to the element, Pr(K) is the polynomial space of degree less than or equal to r.  We apply RT0 space in our coming numerical examples. So our MFE method for approaching the weak formulation for the flow equation is to find ph(·,·)∈Vh and uh(·,·)∈WN0=W0∩Wh+Eh(uB) such that(3.5)(K-1uh,w)-(∇⋅w,ph)=-∫ΓpBw⋅nds,∀w∈Wh,t∈(0,T],(∇⋅uh,v)=(0,v),∀v∈Vh,t∈(0,T].
Because the MFE formulation leads to a saddle point problem for the elliptic equations, we need to use the Mixed-Hybrid algorithms [32] for the pressure equation. Those algorithms mainly use adding unknowns denoting the edge pressure averages, such that the reduced linear system we will solve contains a positive and symmetric definite matrix, and therefore we take the advantages in iterative linear solvers.Based on theRT0 space, (3.5) can be expressed as an algebraic linear system with the unknowns PE (the pressure edge averages)(3.6)ATPE=J.
The derivation of (3.6) is presented in Appendix B.


### 3.4. FVM and FME for Reactive Transport System

#### 3.4.1. Weak Formulation

The weak formulation for the reactive transport system is to find the concentration solutionc∈Vand the diffusive flux solution w∈W0 such that(3.7)(∂(ϕecf)∂t,v)+(∇⋅w,v)+∑K∫∂Kvucf*⋅nds-∑K(cf*u,∇v)=0,∀v∈V,t∈(0,T],(-D-1w,ŵ)+(cf,∇⋅w)=0,∀w∈W0,t∈(0,T],(cf,v)=(cf,0,v),∀v∈V,t=0,
where cf* represents the upwind value of the concentration on an edge. There is the derivation of (3.7) in Appendix C.


#### 3.4.2. MFE Scheme

On the basis of the weak formulation (3.7), the continuous-in-time MFE method to approach the reactive transport system is to find cfh∈Vh and wh∈Wh0 such that(3.8)(∂(ϕecfh)∂t,v)+(∇⋅wh,v)+∑K∫∂Kvucf*⋅nds-∑K(cfh*u,∇v)=0,∀v∈Vh,t∈(0,T],(-D-1wh,ŵ)+(cfh,∇⋅w)=0,∀ŵ∈Wh0,t∈(0,T],(cf,v)=(cf0,v),v∈Vh,t=0.
Then we specify a fully discretized algorithm for the transport. We discretize the simulation time [0,T] into n subintervals:0=t0<t1<⋯<tn-1<tn=T. Let Δtl=tl-tl-1,Δt=maxlΔtl. Provided that there exists a constant C>0 such that Δt≤CminlΔtl, and the transport equation can be solved by semi-implicit Euler method in time together with the combined FVM-MFE method in space.Therefore the fully discretized approximation is to findcfh,l∈Vh, wh,k∈Wh0  (l=0,1,…,n) such that(3.9)(ϕe(cfh,l-cfh,l-1)Δtl,v)+(∇⋅wh,l,v)+∑K∫∂Kvucfh,l-1*⋅nds-∑K(cfh,l-1*u,∇v)=0,∀v∈Vh,(-D-1wh,l,ŵ)+(cfh,l,∇⋅ŵ)=0,∀w∈Wh0,(cfh,0,v)=(cf0,v),∀v∈Vh.
We employ the standard MFE algorithm to solve reactive transport equation (3.7). We here also use RT0 on the reactive transport system, and thus the linear system (2.4) can be specified as the following algebraic linear ordinary differential equation system(3.10)LdXdt+(R⌢+R̂)X=b.
Ordinary differential equations (3.10) are solved by backward Euler method in this paper. The derivation of (3.10) is presented in Appendix D.


## 3.1. Triangular Mesh Generation for Vuggy Fractured Porous Media

Similar to any other numerical modeling, the first step of our algorithm is creating a mesh. The vugs and fractures may be distributed randomly in the domain, so a triangular mesh is required to fit the vuggy fractured porous media much better than a rectangular mesh (though it is easier to generate). A large number of triangular mesh programs are available. Three of them (including DISTMESH, MESHGEN, and TRIANGLE) have much higher mesh quality than any other one in generating adaptive triangular mesh [31]. In this paper, we apply the MESHGEN [32] to generate an adaptive triangular mesh. The fractures are presented in rectangles distributed randomly in the domain considered, and the vugs are described in ellipses distributed randomly.


## 3.2. Mixed Finite Element Method

We employ a mixed finite element (MFE) method for the treatment of the flow equation. MFE method [33] is based on the variational principle that expresses an equilibrium that can be satisfied locally on each finite element. The MFE formulation for the flow equation contains solution for both the scalar variable (pressure) and flux vector (total velocity). Choosing MFE method to approach spaces can satisfy three important properties: local mass conservation, flux continuity, and the same order of convergence (and in some cases super convergence) for both the scalar variable and the flux [33]. MFE method can directly accommodate full permeability tensors and it is more accurate in flux calculation than the conventional finite volume and finite element methods.We first give some symbol notations.  Let(·,·) represent the  L2(Ω) inner product over a domain Ω∈R2 for scalar functions, and (L2(Ω))2 denotes inner product for vector functions. Some necessary spaces are given by(3.1)V=L2(Ω),W=H(div⁡;Ω)={w∈(L2(Ω))2:∇⋅w∈L2(Ω)},W0={w∈H(div⁡;Ω):w⋅n=0on∂Ω},WN0={w∈H(div⁡;Ω):w⋅n=0onΓN}.Based on these spaces, we give the weak formulations for flow equation and reactive equation as follows, respectively.


## 3.3. MFE for Flow Equation

### 3.3.1. Weak Formulation

For the flow equation, the weak formulation is to findp∈Vand u∈WN0+E(uB) such that (3.2)(K-1u,w)-(∇⋅w,p)=-∫ΓpBw⋅nds,∀w∈WN0,t∈(0,T],(∇⋅u,v)=(q,v),∀v∈V,t∈(0,T],
where E(uB) denotes the velocity extension which lets its component agrees with uB on ΓN. There is the derivation of (3.2) in Appendix A.


### 3.3.2. MFE Scheme

We apply the RT (Raviart-Thomas) finite element space [34] to approximate the Darcy velocity. As is well known, the rth order RT space for the two-dimensional triangular element is given by(3.3)Vh(K)=Pr(K),Wh(K)=(x1,Pr(K)x2Pr(K))⊕(Pr(K))2,
where ⊕ represents the direct sum.  We employ the case r=0 in this paper. Therefore, Wh(K) has the form(3.4)Wh(K)={w:w=(aK+bKx1,cK+bKx2),aK,bK,cK∈R}.
Restricted to the element, Pr(K) is the polynomial space of degree less than or equal to r.  We apply RT0 space in our coming numerical examples. So our MFE method for approaching the weak formulation for the flow equation is to find ph(·,·)∈Vh and uh(·,·)∈WN0=W0∩Wh+Eh(uB) such that(3.5)(K-1uh,w)-(∇⋅w,ph)=-∫ΓpBw⋅nds,∀w∈Wh,t∈(0,T],(∇⋅uh,v)=(0,v),∀v∈Vh,t∈(0,T].
Because the MFE formulation leads to a saddle point problem for the elliptic equations, we need to use the Mixed-Hybrid algorithms [32] for the pressure equation. Those algorithms mainly use adding unknowns denoting the edge pressure averages, such that the reduced linear system we will solve contains a positive and symmetric definite matrix, and therefore we take the advantages in iterative linear solvers.Based on theRT0 space, (3.5) can be expressed as an algebraic linear system with the unknowns PE (the pressure edge averages)(3.6)ATPE=J.
The derivation of (3.6) is presented in Appendix B.


## 3.3.1. Weak Formulation

For the flow equation, the weak formulation is to findp∈Vand u∈WN0+E(uB) such that (3.2)(K-1u,w)-(∇⋅w,p)=-∫ΓpBw⋅nds,∀w∈WN0,t∈(0,T],(∇⋅u,v)=(q,v),∀v∈V,t∈(0,T],
where E(uB) denotes the velocity extension which lets its component agrees with uB on ΓN. There is the derivation of (3.2) in Appendix A.


## 3.3.2. MFE Scheme

We apply the RT (Raviart-Thomas) finite element space [34] to approximate the Darcy velocity. As is well known, the rth order RT space for the two-dimensional triangular element is given by(3.3)Vh(K)=Pr(K),Wh(K)=(x1,Pr(K)x2Pr(K))⊕(Pr(K))2,
where ⊕ represents the direct sum.  We employ the case r=0 in this paper. Therefore, Wh(K) has the form(3.4)Wh(K)={w:w=(aK+bKx1,cK+bKx2),aK,bK,cK∈R}.
Restricted to the element, Pr(K) is the polynomial space of degree less than or equal to r.  We apply RT0 space in our coming numerical examples. So our MFE method for approaching the weak formulation for the flow equation is to find ph(·,·)∈Vh and uh(·,·)∈WN0=W0∩Wh+Eh(uB) such that(3.5)(K-1uh,w)-(∇⋅w,ph)=-∫ΓpBw⋅nds,∀w∈Wh,t∈(0,T],(∇⋅uh,v)=(0,v),∀v∈Vh,t∈(0,T].
Because the MFE formulation leads to a saddle point problem for the elliptic equations, we need to use the Mixed-Hybrid algorithms [32] for the pressure equation. Those algorithms mainly use adding unknowns denoting the edge pressure averages, such that the reduced linear system we will solve contains a positive and symmetric definite matrix, and therefore we take the advantages in iterative linear solvers.Based on theRT0 space, (3.5) can be expressed as an algebraic linear system with the unknowns PE (the pressure edge averages)(3.6)ATPE=J.
The derivation of (3.6) is presented in Appendix B.


## 3.4. FVM and FME for Reactive Transport System

### 3.4.1. Weak Formulation

The weak formulation for the reactive transport system is to find the concentration solutionc∈Vand the diffusive flux solution w∈W0 such that(3.7)(∂(ϕecf)∂t,v)+(∇⋅w,v)+∑K∫∂Kvucf*⋅nds-∑K(cf*u,∇v)=0,∀v∈V,t∈(0,T],(-D-1w,ŵ)+(cf,∇⋅w)=0,∀w∈W0,t∈(0,T],(cf,v)=(cf,0,v),∀v∈V,t=0,
where cf* represents the upwind value of the concentration on an edge. There is the derivation of (3.7) in Appendix C.


### 3.4.2. MFE Scheme

On the basis of the weak formulation (3.7), the continuous-in-time MFE method to approach the reactive transport system is to find cfh∈Vh and wh∈Wh0 such that(3.8)(∂(ϕecfh)∂t,v)+(∇⋅wh,v)+∑K∫∂Kvucf*⋅nds-∑K(cfh*u,∇v)=0,∀v∈Vh,t∈(0,T],(-D-1wh,ŵ)+(cfh,∇⋅w)=0,∀ŵ∈Wh0,t∈(0,T],(cf,v)=(cf0,v),v∈Vh,t=0.
Then we specify a fully discretized algorithm for the transport. We discretize the simulation time [0,T] into n subintervals:0=t0<t1<⋯<tn-1<tn=T. Let Δtl=tl-tl-1,Δt=maxlΔtl. Provided that there exists a constant C>0 such that Δt≤CminlΔtl, and the transport equation can be solved by semi-implicit Euler method in time together with the combined FVM-MFE method in space.Therefore the fully discretized approximation is to findcfh,l∈Vh, wh,k∈Wh0  (l=0,1,…,n) such that(3.9)(ϕe(cfh,l-cfh,l-1)Δtl,v)+(∇⋅wh,l,v)+∑K∫∂Kvucfh,l-1*⋅nds-∑K(cfh,l-1*u,∇v)=0,∀v∈Vh,(-D-1wh,l,ŵ)+(cfh,l,∇⋅ŵ)=0,∀w∈Wh0,(cfh,0,v)=(cf0,v),∀v∈Vh.
We employ the standard MFE algorithm to solve reactive transport equation (3.7). We here also use RT0 on the reactive transport system, and thus the linear system (2.4) can be specified as the following algebraic linear ordinary differential equation system(3.10)LdXdt+(R⌢+R̂)X=b.
Ordinary differential equations (3.10) are solved by backward Euler method in this paper. The derivation of (3.10) is presented in Appendix D.


## 3.4.1. Weak Formulation

The weak formulation for the reactive transport system is to find the concentration solutionc∈Vand the diffusive flux solution w∈W0 such that(3.7)(∂(ϕecf)∂t,v)+(∇⋅w,v)+∑K∫∂Kvucf*⋅nds-∑K(cf*u,∇v)=0,∀v∈V,t∈(0,T],(-D-1w,ŵ)+(cf,∇⋅w)=0,∀w∈W0,t∈(0,T],(cf,v)=(cf,0,v),∀v∈V,t=0,
where cf* represents the upwind value of the concentration on an edge. There is the derivation of (3.7) in Appendix C.


## 3.4.2. MFE Scheme

On the basis of the weak formulation (3.7), the continuous-in-time MFE method to approach the reactive transport system is to find cfh∈Vh and wh∈Wh0 such that(3.8)(∂(ϕecfh)∂t,v)+(∇⋅wh,v)+∑K∫∂Kvucf*⋅nds-∑K(cfh*u,∇v)=0,∀v∈Vh,t∈(0,T],(-D-1wh,ŵ)+(cfh,∇⋅w)=0,∀ŵ∈Wh0,t∈(0,T],(cf,v)=(cf0,v),v∈Vh,t=0.
Then we specify a fully discretized algorithm for the transport. We discretize the simulation time [0,T] into n subintervals:0=t0<t1<⋯<tn-1<tn=T. Let Δtl=tl-tl-1,Δt=maxlΔtl. Provided that there exists a constant C>0 such that Δt≤CminlΔtl, and the transport equation can be solved by semi-implicit Euler method in time together with the combined FVM-MFE method in space.Therefore the fully discretized approximation is to findcfh,l∈Vh, wh,k∈Wh0  (l=0,1,…,n) such that(3.9)(ϕe(cfh,l-cfh,l-1)Δtl,v)+(∇⋅wh,l,v)+∑K∫∂Kvucfh,l-1*⋅nds-∑K(cfh,l-1*u,∇v)=0,∀v∈Vh,(-D-1wh,l,ŵ)+(cfh,l,∇⋅ŵ)=0,∀w∈Wh0,(cfh,0,v)=(cf0,v),∀v∈Vh.
We employ the standard MFE algorithm to solve reactive transport equation (3.7). We here also use RT0 on the reactive transport system, and thus the linear system (2.4) can be specified as the following algebraic linear ordinary differential equation system(3.10)LdXdt+(R⌢+R̂)X=b.
Ordinary differential equations (3.10) are solved by backward Euler method in this paper. The derivation of (3.10) is presented in Appendix D.


## 4. Numerical Results

### 4.1. Simulation Cases

We mainly simulate two sequential phases for an infiltration process with different initial and boundary conditions: at first, the clean water infiltration process, with initial concentrationc1(x,0)=0 and boundary condition in the inflow boundary c1(x,t)=0, x∈Γinflow. In the first phase, the polluted water infiltration process, with initial concentration in the whole region c2(x,0)=0.0 and boundary condition in the inflow boundaryc1(x,t)=1.0  x∈Γinflow. In the second phase, again the clean water infiltration process, with initial concentration in the whole region c3(x,0)=c1(x,T1) and boundary condition in the inflow boundary c3(x,t)=0, x∈Γinflow, where c1(x,T1) is the concentration distribution at T1 in the first phase. The domain considered for all cases in this paper is a bounded rectangular domain Ω⊂R2 of (0,0.80m)×(0,0.80m) with randomly generated vugs and fractures into an adaptive triangular mesh. The mesh is fixed for all time during the simulation. The region around the fractures and/or vugs is deeply refined because of dramatically changing of some parameters for the models. For all cases in the paper, a series of standard parameters applied in practice are listed in Table 1 as follows [1].Table 1
Standard parameters for the mathematical model in practice.
Saturated hydraulic conductivity,Ks (cm/s)Saturated effective diffusion coefficient,De (cm2/s)Effective porosityDistribution coefficientKd1.0e−125.0e−110.1840.395Case 0 as a base case describes the porous media without vugs and fractures to be compared with the other cases, that is, the fractured vuggy porous media with vugs and fractures distributed differently. In Case 1, in the fractured vuggy porous media, the vugs are located on the inflow boundary connected with the fractures, and the fracture does not elongate to the outflow boundary. In Case 2, the vugs are still located on the inflow boundary connected with the fracture, but the fractures extend to the outflow boundary. In Case 3, the vugs connected with fracture are located within the region, and both the vugs and the fractures do not elongate to any boundary. In Cases 5 and 6, the vugs lie on the outflow boundary connected with fracture, but in Case 5 the fractures are connected with the inflow boundary and the fractures are not in Case 6. The vugs in Case 7 are connected with fractures, and the fractures are connected with inflow boundary and outflow boundary. In the last case, the fractures are located within the region connecting the vugs that lie on the inflow boundary and the outflow boundary. The no-flow boundaries in all cases are taken as the top (y=0.80 m)  and the bottom boundaries (y=0.80 m). The inflow boundary and outflow boundary are taken as the left boundary (x=0.0 m) and the right boundary (x=0.80 m) of the domain, respectively.  A constant pressure of 0 (the gauge pressure against a reference pressure) is specified on the outflow boundary (x=0.80 m). In the first phase of the simulation, the polluted water is continuously injected on the left boundary (x=0 m), that is, the inflow boundary, and the region is initially full of clean water. In the second phase of the simulation, the clean water is continuously injected at the final of the first phase. We simulate up to 35,000 years and 35,000 for the first phase and the second phase for all cases, respectively. The injection on the inflow boundary for the two phases is showed in Figure 2.Figure 2
Injection on the inflow boundary for the two phases of the simulation.


### 4.2. Simulation Results

Case 0: Matrix Porous Media
In this situation, the area is matrix porous media as a base case to be compared with 8 other cases. In the first phase the polluted water is constantly injected into the inflow boundary at the time interval(0,T1]. The conductivity distribution is showed in Figure 3. An adaptive triangular mesh for this domain is generated (Figure 3).  We approximate the Darcy equation and transport equation using RT0-MFE and semi-implicit FVM-MFE presented before, respectively, with a uniform time step of 350 years. The pressure field and the streamlines field (Figure 4) are obtained. Figure 5 demonstrates the results of simulated profiles at different time. The results show that the polluted water is uniformly diffusing from the inflow boundary through the domain to the outflow boundary in the first phase. In the second phase, the clean water is being injected on the inflow boundary. Figure 6 shows the results of the simulated profiles at different time for the second phase. The results demonstrate that the polluted water is uniformly ejected out from the outflow boundary. The effluent concentration on the outlet is obtained for the two phases to be compared with the other cases.Conductivity distribution and adaptive triangular mesh for Case 0.
(a)
Conductivity distribution(b)
Adaptive triangular meshPressure distribution and streamlines field for base case.
(a)
Pressure distribution(b)
Streamlines fieldConcentration distribution at different time for the first phase for base case.
(a)
The 3500th year(b)
The 14000th year(c)
The 35000th yearConcentration distribution at different time for the second phase for base case.
(a)
The 12000 year(b)
The 15000 year(c)
The 20000 yearCase 1: Vugs on the Inflow Boundary Connected with Fractures inside the Domain
In this case, vugs are located on the inflow boundary of the considered domain connected with fractures not touching the outflow boundary, which is depicted in Figure7. The pressure distribution and the streamlines field are shown in Figure 8.  From Figure 8, it is obviously demonstrated that the vugs and fractures affect the pressure distribution compared with the base case. The vugs and fractures have corresponding locally irregular pressure, and the pressure around vugs and fractures is higher than other regions in the domain. Therefore the quantity of velocity in the vugs and fractures is much larger than that in the matrix media. Furthermore, we can observe that the streamlines are inclined to converging into the vugs and the fractures near the inflow part but diverge from the fractures on the right part, which show that the vugs and fractures are the main pathways for contaminant transport by convection. Results of simulated concentration sections are shown in Figure 9 at different time. The results show that the contaminant transports essentially via convection in the domain of vugs and fractures from the first year to the 3500th year. After the 14000th year, convection and diffusion through the matrix area also begin to take obvious effect on the contaminant transport behavior. During the second phase of the simulation, from the 35001th to 38500th year the vuggy and fractured region is first getting much cleaner than that in the matrix, which is showed in Figure 10.  After the 49000th year, the clean water begins to infiltrate into the matrix porous media displayed in Figure 10.Conductivity distribution and adaptive triangular mesh for Case 1.
(a)
Conductivity distribution(b)
Adaptive triangular meshPressure distribution, streamlines field, and average effluent concentration for Case 1.
(a)
Pressure distribution(b)
Streamlines field(c)
Average effluent concentrationConcentration distribution at different time for the first phase for Case 1.
(a)
The 3500th year(b)
The 14000 year(c)
The 35000th yearConcentration distribution at different time for the second phase for Case 1.
(a)
The 38500th year(b)
The 49000th year(c)
The 70000th yearCase 2: Vugs on the Inflow Boundary Connected with Fractures Extending to the Outflow Boundary
Compared with Case 1, the fractures in this case are extending to the outflow boundary, the conductivity and the adaptive triangular mesh of which are shown in Figure11. The pressure distribution and the streamlines field are shown in Figure 12. Similar to Case 1, the vugs and fractures affect the pressure distribution compared with the base case, and the vugs and fractures have corresponding locally irregular pressure.  But different from Case 1, the streamlines are inclined to converging into the vugs and the fractures from the inlet to the outlet, which show that the vugs and fractures are still the main pathways for transporting contaminant by convection. Results of simulated concentration profiles are shown in Figure 13 at different time for the first phase. From the first year to the 14000th year, the contaminant transports mainly through the vugs and the fractures because of the much faster velocity in the vugs and the fractures than that in matrix. And the pollutant quickly passes through fractures to the outflow boundary. From the 14001th year to the 35000th year, the matrix begins to become a significant role in contaminant transport, as is shown in Figure 13. During the second phase of the simulation, the clean water quickly passes through the vugs and fractures to the outlet from the 35001th year to the 38500th year demonstrated in Figure 14, and the matrix near the inflow boundary begins getting clean from the 38501th year, which is shown in Figure 14.Conductivity distribution and adaptive triangular mesh for Case 2.
(a)
Conductivity distribution(b)
Adaptive triangular meshPressure distribution, streamlines field, and average effluent concentration for Case 2.
(a)
Pressure distribution(b)
Streamlines field(c)
Average effluent concentrationConcentration distribution at different time for the first phase for Case 2.
(a)
The 3500th year(b)
The 14000th year(c)
The 35000th yearConcentration distribution at different time for the second phase for Case 2.
(a)
The 38500th year(b)
The 49000th year(c)
The 70000th yearCase 3: Vugs Connected with Fractures inside the Domain
In this case, none of the vugs and fractures touches the boundary, the conductivity and adaptive triangular mesh of which are displayed in Figure15. We can observe that the vugs observably affect the pressure distribution from the pressure field in Figure 16. On the contrary, the fractures do not affect the pressure distribution so obviously as the vugs. From the streamlines field in Figure 16, we also find out that the streamlines converge into the vugs and fractures for the inflow boundary and diverge from the endpoint of the fractures. For the first phase of the simulation, the results of simulated concentration profiles are shown in Figure 17. From the initial situation to 1400th year, the contaminant mainly transports in the matrix porous media via convection and diffusion. From the 1401th year to the 3500th year, the contaminant mainly infiltrates into the vugs and fractures (Figure 17). From the 3501th year, the matrix starts to significantly take effect on the transport through diffusion and convection. When injecting clean water in the second phase, from the 35001th year the porous matrix media are firstly getting cleaner than that in matrix, and from the 38500th year the clean water quickly infiltrate into the vugs and fractures with the contaminant remained in the matrix media. From 49000th year the clean water starts to infiltrate the matrix (Figure 18).Conductivity distribution and adaptive triangular mesh for Case 3.
(a)
Conductivity distribution(b)
Adaptive triangular meshPressure distribution, streamlines field, and average effluent concentration for Case 3.
(a)
Pressure distribution(b)
Streamlines field(c)
Average effluent concentrationConcentration distribution at different time for the first phase for Case 3.
(a)
The 3500th year(b)
The 14000th year(c)
The 35000th yearConcentration distribution at different time for the second phase for Case 3.
(a)
The 38500th year(b)
The 49000th year(c)
The 70000th yearCase 4: Vugs on the Outflow Boundary Connected with Fractures inside the Domain
Now we consider the vugs on the outlet connected with fractures no touching the inlet. The conductivity distribution and adaptive triangular mesh of the case are displayed in Figure19. It is observed that the pressure field (Figure 20) is not affected by vugs so significantly as that in Case 3. This is because we take the pressure a 0 on the outflow boundary. From the streamlines field (Figure 20), we can find out that the streamlines converge into the fractures from the inflow boundary and diverge at the end of the fractures. From the results of the simulated concentration sections (Figure 21). We can observe that from first year to the 3500th year, the contaminant transports mainly through the matrix via diffusion and convection. From the 3501th year to the 5000th year, the fractures and vugs begin to influence the transport of contaminant via convection (Figure 21). From the 5001th year to the 35000th year the contaminant transports through the whole fractured vuggy porous media (Figure 21). In the second phase of the simulation, the clean water almost uniformly flows through the matrix from the 35001th year to the 38500th year (Figure 22), and from the 3900th year to the 44000th year, the clean water mainly passes through the vugs and fractures (Figure 22). From the 44001th year to 70000th year, clean water begins to pass the whole fractured vuggy porous media (Figure 22).Conductivity distribution and adaptive triangular mesh for Case 4.
(a)
Conductivity distribution(b)
Adaptive triangular meshPressure distribution, streamlines field, and average effluent concentration for Case 4.
(a)
Pressure distribution(b)
Streamlines field(c)
Average effluent concentrationConcentration distribution at different time for the first phase for Case 4.
(a)
The 3500th year(b)
The 14000th year(c)
The 35000th yearConcentration distribution at different time for the second phase for Case 4.
(a)
The 38500th year(b)
The 49000th year(c)
The 70000th yearCase 5:Vugs on the Outflow Boundary Connected with Fractures Extending to the Inflow Boundary
Similar to Case 4, but the fractures are now extending to the inflow boundary. The conductivity distribution and adaptive triangular mesh of this situation are given in Figure23. The pressure and streamlines fields are displayed in Figure 24. The pressure field shows that the fractures and vugs do not affect the pressure distribution. But the fractures significantly affect the streamlines distribution (Figure 24). From the results of the simulated concentration, different from Case 4, the contaminant transport mainly passes through vugs and fractures via convection from year 1 to year 3500 in the first phase, and the contaminant quickly passes through vugs and fractures to the outlet (Figure 25). From the 3501th year the contaminant starts to infiltrate into the matrix media (Figure 25). In the second phase, the clean water firstly passes through the fractures and vugs (Figure 26).Conductivity distribution and adaptive triangular mesh for Case 5.
(a)
Conductivity distribution(b)
Adaptive triangular meshPressure distribution, streamlines field, and average effluent concentration for Case 5.
(a)
Pressure distribution(b)
Streamlines field(c)
Average effluent concentrationConcentration distribution at different time for the first phase for Case 5.
(a)
The 3500th year(b)
The 14000th year(c)
The 35000th yearConcentration distribution at different time for the second phase for Case 5.
(a)
The 38500th year(b)
The 49000th year(c)
The 70000 yearCase 6: Vugs inside the Domain Connected with Fractures Extending Both Inflow and Outflow Boundaries
Now the vugs are located inside the domain connected with fractures extending to the inlet and outlet. The conductivity distribution and adaptive triangular mesh are displayed in Figure27. The pressure field and streamlines field are demonstrated in Figure 28. The pressure field suggests that the fractures near the inflow boundary and vugs significantly affect the pressure distribution while the fractures near the outflow boundary not. The streamlines field shows that the streamlines converge into the fractures and the vugs. The results of the simulated concentration are showed in Figure 29. From Figure 29, we can observe that the contaminant quickly passes through the vugs and fractures to the outlet via convection during year 1 to year 3500, after that the contaminant starts to pass through the matrix. In the second phase, the clean water behaves as the contaminant in the first phase as shown in Figure 30.Conductivity distribution and adaptive triangular mesh for Case 6.
(a)
Conductivity distribution(b)
Adaptive triangular meshPressure distribution, streamlines field, and average effluent concentration for Case 6.
(a)
Pressure distribution(b)
Streamlines field(c)
Average effluent concentrationConcentration distribution at different time for the first phase for Case 6.
(a)
The 3500th year(b)
The 14000th year(c)
The 35000th yearConcentration distribution at different time for the second phase for Case 6.
(a)
The 38500th year(b)
The 49000th year(c)
The 70000th yearCase 7: Vugs on the Inflow Boundary and Outflow Boundary Connected with Fractures through the Domain
In this case, the vugs are located on the inlet and outlet connected with fractures, the conductivity distribution and adaptive triangular mesh of which are displayed in Figure31. The pressure and streamlines fields are showed in Figure 32. From Figure 32 we can observe that the vugs near the inflow boundary obviously affect the pressure distribution while the fractures and the vugs near the outflow boundary not. From year 1 to year 3500, the contaminant transports mainly through the vugs and fractures via convection (Figure 33). From the 3501th year, the contaminant starts to infiltrate into the matrix porous media (Figure 33). In the second phase, the clean water quickly passes through the vugs and fractures to the outlet during year 35001 and the 38500th year. After that the clean water starts to infiltrate into the whole domain (Figure 34).Conductivity distribution and adaptive triangular mesh for Case 7.
(a)
Conductivity distribution(b)
Adaptive triangular meshPressure distribution, streamlines field, and Average effluent concentration for Case 7.
(a)
Pressure distribution(b)
Streamlines field(c)
Average effluent concentrationConcentration distribution at different time for the first phase for Case 7.
(a)
The 3500th year(b)
The 14000th year(c)
The 35000th yearConcentration distribution at different time for the second phase for Case 7.
(a)
The 38500th year(b)
The 49000th year(c)
The 70000 year


### 4.3. Analysis of the Concentration on the Outflow Boundary at Different Time

We get the effluent concentration on the outlet for every case at different time (Figures8, 12, 16, 20, 24, 28, and 32). For Case 1 to Case 7, they have a common feature that the effluent concentration quickly becomes higher in shorter time than that in the base case during the first phase, and then the concentration increases gradually until it gets to the highest (1.0). This suggests that there are time multiscale phenomena in the contaminant transport in karstic aquifers.  A main reason may be that the vugs and fractures have much higher porosity than that in matrix porous media, and the velocity in the vugs and fractures is also higher than that in matrix porous media. Once the clean water is injected into the inlet, the fractured vuggy porous media, the effluent concentration becomes lower in shorter time than that in the base case except for the Case 4. It is easily to understand that the vugs and fractures are the main pathways for contaminant transport. But it takes longer time for the vuggy fractured porous media to get the lowest concentration than that in base case (Figures 8, 12, 16, 20, 24, 28, and 32). It may be because when the clean water is injected into the inlet, it firstly infiltrates into the vugs and fractures, but there remains much contaminant in the matrix porous media by adsorption. So we can conclude that karst groundwater becomes polluted more easily and in shorter time periods than that in non-karstic aquifers. Even though the aquifers are injected clean water which has been polluted by some contaminants, it will take a much longer time to render the groundwater to be clean again than that in non-karstic areas, because of the significant difference of adsorption between the vugs, fractures, and matrix porous media.


## 4.1. Simulation Cases

We mainly simulate two sequential phases for an infiltration process with different initial and boundary conditions: at first, the clean water infiltration process, with initial concentrationc1(x,0)=0 and boundary condition in the inflow boundary c1(x,t)=0, x∈Γinflow. In the first phase, the polluted water infiltration process, with initial concentration in the whole region c2(x,0)=0.0 and boundary condition in the inflow boundaryc1(x,t)=1.0  x∈Γinflow. In the second phase, again the clean water infiltration process, with initial concentration in the whole region c3(x,0)=c1(x,T1) and boundary condition in the inflow boundary c3(x,t)=0, x∈Γinflow, where c1(x,T1) is the concentration distribution at T1 in the first phase. The domain considered for all cases in this paper is a bounded rectangular domain Ω⊂R2 of (0,0.80m)×(0,0.80m) with randomly generated vugs and fractures into an adaptive triangular mesh. The mesh is fixed for all time during the simulation. The region around the fractures and/or vugs is deeply refined because of dramatically changing of some parameters for the models. For all cases in the paper, a series of standard parameters applied in practice are listed in Table 1 as follows [1].Table 1
Standard parameters for the mathematical model in practice.
Saturated hydraulic conductivity,Ks (cm/s)Saturated effective diffusion coefficient,De (cm2/s)Effective porosityDistribution coefficientKd1.0e−125.0e−110.1840.395Case 0 as a base case describes the porous media without vugs and fractures to be compared with the other cases, that is, the fractured vuggy porous media with vugs and fractures distributed differently. In Case 1, in the fractured vuggy porous media, the vugs are located on the inflow boundary connected with the fractures, and the fracture does not elongate to the outflow boundary. In Case 2, the vugs are still located on the inflow boundary connected with the fracture, but the fractures extend to the outflow boundary. In Case 3, the vugs connected with fracture are located within the region, and both the vugs and the fractures do not elongate to any boundary. In Cases 5 and 6, the vugs lie on the outflow boundary connected with fracture, but in Case 5 the fractures are connected with the inflow boundary and the fractures are not in Case 6. The vugs in Case 7 are connected with fractures, and the fractures are connected with inflow boundary and outflow boundary. In the last case, the fractures are located within the region connecting the vugs that lie on the inflow boundary and the outflow boundary. The no-flow boundaries in all cases are taken as the top (y=0.80 m)  and the bottom boundaries (y=0.80 m). The inflow boundary and outflow boundary are taken as the left boundary (x=0.0 m) and the right boundary (x=0.80 m) of the domain, respectively.  A constant pressure of 0 (the gauge pressure against a reference pressure) is specified on the outflow boundary (x=0.80 m). In the first phase of the simulation, the polluted water is continuously injected on the left boundary (x=0 m), that is, the inflow boundary, and the region is initially full of clean water. In the second phase of the simulation, the clean water is continuously injected at the final of the first phase. We simulate up to 35,000 years and 35,000 for the first phase and the second phase for all cases, respectively. The injection on the inflow boundary for the two phases is showed in Figure 2.Figure 2
Injection on the inflow boundary for the two phases of the simulation.


## 4.2. Simulation Results

Case 0: Matrix Porous Media
In this situation, the area is matrix porous media as a base case to be compared with 8 other cases. In the first phase the polluted water is constantly injected into the inflow boundary at the time interval(0,T1]. The conductivity distribution is showed in Figure 3. An adaptive triangular mesh for this domain is generated (Figure 3).  We approximate the Darcy equation and transport equation using RT0-MFE and semi-implicit FVM-MFE presented before, respectively, with a uniform time step of 350 years. The pressure field and the streamlines field (Figure 4) are obtained. Figure 5 demonstrates the results of simulated profiles at different time. The results show that the polluted water is uniformly diffusing from the inflow boundary through the domain to the outflow boundary in the first phase. In the second phase, the clean water is being injected on the inflow boundary. Figure 6 shows the results of the simulated profiles at different time for the second phase. The results demonstrate that the polluted water is uniformly ejected out from the outflow boundary. The effluent concentration on the outlet is obtained for the two phases to be compared with the other cases.Conductivity distribution and adaptive triangular mesh for Case 0.
(a)
Conductivity distribution(b)
Adaptive triangular meshPressure distribution and streamlines field for base case.
(a)
Pressure distribution(b)
Streamlines fieldConcentration distribution at different time for the first phase for base case.
(a)
The 3500th year(b)
The 14000th year(c)
The 35000th yearConcentration distribution at different time for the second phase for base case.
(a)
The 12000 year(b)
The 15000 year(c)
The 20000 yearCase 1: Vugs on the Inflow Boundary Connected with Fractures inside the Domain
In this case, vugs are located on the inflow boundary of the considered domain connected with fractures not touching the outflow boundary, which is depicted in Figure7. The pressure distribution and the streamlines field are shown in Figure 8.  From Figure 8, it is obviously demonstrated that the vugs and fractures affect the pressure distribution compared with the base case. The vugs and fractures have corresponding locally irregular pressure, and the pressure around vugs and fractures is higher than other regions in the domain. Therefore the quantity of velocity in the vugs and fractures is much larger than that in the matrix media. Furthermore, we can observe that the streamlines are inclined to converging into the vugs and the fractures near the inflow part but diverge from the fractures on the right part, which show that the vugs and fractures are the main pathways for contaminant transport by convection. Results of simulated concentration sections are shown in Figure 9 at different time. The results show that the contaminant transports essentially via convection in the domain of vugs and fractures from the first year to the 3500th year. After the 14000th year, convection and diffusion through the matrix area also begin to take obvious effect on the contaminant transport behavior. During the second phase of the simulation, from the 35001th to 38500th year the vuggy and fractured region is first getting much cleaner than that in the matrix, which is showed in Figure 10.  After the 49000th year, the clean water begins to infiltrate into the matrix porous media displayed in Figure 10.Conductivity distribution and adaptive triangular mesh for Case 1.
(a)
Conductivity distribution(b)
Adaptive triangular meshPressure distribution, streamlines field, and average effluent concentration for Case 1.
(a)
Pressure distribution(b)
Streamlines field(c)
Average effluent concentrationConcentration distribution at different time for the first phase for Case 1.
(a)
The 3500th year(b)
The 14000 year(c)
The 35000th yearConcentration distribution at different time for the second phase for Case 1.
(a)
The 38500th year(b)
The 49000th year(c)
The 70000th yearCase 2: Vugs on the Inflow Boundary Connected with Fractures Extending to the Outflow Boundary
Compared with Case 1, the fractures in this case are extending to the outflow boundary, the conductivity and the adaptive triangular mesh of which are shown in Figure11. The pressure distribution and the streamlines field are shown in Figure 12. Similar to Case 1, the vugs and fractures affect the pressure distribution compared with the base case, and the vugs and fractures have corresponding locally irregular pressure.  But different from Case 1, the streamlines are inclined to converging into the vugs and the fractures from the inlet to the outlet, which show that the vugs and fractures are still the main pathways for transporting contaminant by convection. Results of simulated concentration profiles are shown in Figure 13 at different time for the first phase. From the first year to the 14000th year, the contaminant transports mainly through the vugs and the fractures because of the much faster velocity in the vugs and the fractures than that in matrix. And the pollutant quickly passes through fractures to the outflow boundary. From the 14001th year to the 35000th year, the matrix begins to become a significant role in contaminant transport, as is shown in Figure 13. During the second phase of the simulation, the clean water quickly passes through the vugs and fractures to the outlet from the 35001th year to the 38500th year demonstrated in Figure 14, and the matrix near the inflow boundary begins getting clean from the 38501th year, which is shown in Figure 14.Conductivity distribution and adaptive triangular mesh for Case 2.
(a)
Conductivity distribution(b)
Adaptive triangular meshPressure distribution, streamlines field, and average effluent concentration for Case 2.
(a)
Pressure distribution(b)
Streamlines field(c)
Average effluent concentrationConcentration distribution at different time for the first phase for Case 2.
(a)
The 3500th year(b)
The 14000th year(c)
The 35000th yearConcentration distribution at different time for the second phase for Case 2.
(a)
The 38500th year(b)
The 49000th year(c)
The 70000th yearCase 3: Vugs Connected with Fractures inside the Domain
In this case, none of the vugs and fractures touches the boundary, the conductivity and adaptive triangular mesh of which are displayed in Figure15. We can observe that the vugs observably affect the pressure distribution from the pressure field in Figure 16. On the contrary, the fractures do not affect the pressure distribution so obviously as the vugs. From the streamlines field in Figure 16, we also find out that the streamlines converge into the vugs and fractures for the inflow boundary and diverge from the endpoint of the fractures. For the first phase of the simulation, the results of simulated concentration profiles are shown in Figure 17. From the initial situation to 1400th year, the contaminant mainly transports in the matrix porous media via convection and diffusion. From the 1401th year to the 3500th year, the contaminant mainly infiltrates into the vugs and fractures (Figure 17). From the 3501th year, the matrix starts to significantly take effect on the transport through diffusion and convection. When injecting clean water in the second phase, from the 35001th year the porous matrix media are firstly getting cleaner than that in matrix, and from the 38500th year the clean water quickly infiltrate into the vugs and fractures with the contaminant remained in the matrix media. From 49000th year the clean water starts to infiltrate the matrix (Figure 18).Conductivity distribution and adaptive triangular mesh for Case 3.
(a)
Conductivity distribution(b)
Adaptive triangular meshPressure distribution, streamlines field, and average effluent concentration for Case 3.
(a)
Pressure distribution(b)
Streamlines field(c)
Average effluent concentrationConcentration distribution at different time for the first phase for Case 3.
(a)
The 3500th year(b)
The 14000th year(c)
The 35000th yearConcentration distribution at different time for the second phase for Case 3.
(a)
The 38500th year(b)
The 49000th year(c)
The 70000th yearCase 4: Vugs on the Outflow Boundary Connected with Fractures inside the Domain
Now we consider the vugs on the outlet connected with fractures no touching the inlet. The conductivity distribution and adaptive triangular mesh of the case are displayed in Figure19. It is observed that the pressure field (Figure 20) is not affected by vugs so significantly as that in Case 3. This is because we take the pressure a 0 on the outflow boundary. From the streamlines field (Figure 20), we can find out that the streamlines converge into the fractures from the inflow boundary and diverge at the end of the fractures. From the results of the simulated concentration sections (Figure 21). We can observe that from first year to the 3500th year, the contaminant transports mainly through the matrix via diffusion and convection. From the 3501th year to the 5000th year, the fractures and vugs begin to influence the transport of contaminant via convection (Figure 21). From the 5001th year to the 35000th year the contaminant transports through the whole fractured vuggy porous media (Figure 21). In the second phase of the simulation, the clean water almost uniformly flows through the matrix from the 35001th year to the 38500th year (Figure 22), and from the 3900th year to the 44000th year, the clean water mainly passes through the vugs and fractures (Figure 22). From the 44001th year to 70000th year, clean water begins to pass the whole fractured vuggy porous media (Figure 22).Conductivity distribution and adaptive triangular mesh for Case 4.
(a)
Conductivity distribution(b)
Adaptive triangular meshPressure distribution, streamlines field, and average effluent concentration for Case 4.
(a)
Pressure distribution(b)
Streamlines field(c)
Average effluent concentrationConcentration distribution at different time for the first phase for Case 4.
(a)
The 3500th year(b)
The 14000th year(c)
The 35000th yearConcentration distribution at different time for the second phase for Case 4.
(a)
The 38500th year(b)
The 49000th year(c)
The 70000th yearCase 5:Vugs on the Outflow Boundary Connected with Fractures Extending to the Inflow Boundary
Similar to Case 4, but the fractures are now extending to the inflow boundary. The conductivity distribution and adaptive triangular mesh of this situation are given in Figure23. The pressure and streamlines fields are displayed in Figure 24. The pressure field shows that the fractures and vugs do not affect the pressure distribution. But the fractures significantly affect the streamlines distribution (Figure 24). From the results of the simulated concentration, different from Case 4, the contaminant transport mainly passes through vugs and fractures via convection from year 1 to year 3500 in the first phase, and the contaminant quickly passes through vugs and fractures to the outlet (Figure 25). From the 3501th year the contaminant starts to infiltrate into the matrix media (Figure 25). In the second phase, the clean water firstly passes through the fractures and vugs (Figure 26).Conductivity distribution and adaptive triangular mesh for Case 5.
(a)
Conductivity distribution(b)
Adaptive triangular meshPressure distribution, streamlines field, and average effluent concentration for Case 5.
(a)
Pressure distribution(b)
Streamlines field(c)
Average effluent concentrationConcentration distribution at different time for the first phase for Case 5.
(a)
The 3500th year(b)
The 14000th year(c)
The 35000th yearConcentration distribution at different time for the second phase for Case 5.
(a)
The 38500th year(b)
The 49000th year(c)
The 70000 yearCase 6: Vugs inside the Domain Connected with Fractures Extending Both Inflow and Outflow Boundaries
Now the vugs are located inside the domain connected with fractures extending to the inlet and outlet. The conductivity distribution and adaptive triangular mesh are displayed in Figure27. The pressure field and streamlines field are demonstrated in Figure 28. The pressure field suggests that the fractures near the inflow boundary and vugs significantly affect the pressure distribution while the fractures near the outflow boundary not. The streamlines field shows that the streamlines converge into the fractures and the vugs. The results of the simulated concentration are showed in Figure 29. From Figure 29, we can observe that the contaminant quickly passes through the vugs and fractures to the outlet via convection during year 1 to year 3500, after that the contaminant starts to pass through the matrix. In the second phase, the clean water behaves as the contaminant in the first phase as shown in Figure 30.Conductivity distribution and adaptive triangular mesh for Case 6.
(a)
Conductivity distribution(b)
Adaptive triangular meshPressure distribution, streamlines field, and average effluent concentration for Case 6.
(a)
Pressure distribution(b)
Streamlines field(c)
Average effluent concentrationConcentration distribution at different time for the first phase for Case 6.
(a)
The 3500th year(b)
The 14000th year(c)
The 35000th yearConcentration distribution at different time for the second phase for Case 6.
(a)
The 38500th year(b)
The 49000th year(c)
The 70000th yearCase 7: Vugs on the Inflow Boundary and Outflow Boundary Connected with Fractures through the Domain
In this case, the vugs are located on the inlet and outlet connected with fractures, the conductivity distribution and adaptive triangular mesh of which are displayed in Figure31. The pressure and streamlines fields are showed in Figure 32. From Figure 32 we can observe that the vugs near the inflow boundary obviously affect the pressure distribution while the fractures and the vugs near the outflow boundary not. From year 1 to year 3500, the contaminant transports mainly through the vugs and fractures via convection (Figure 33). From the 3501th year, the contaminant starts to infiltrate into the matrix porous media (Figure 33). In the second phase, the clean water quickly passes through the vugs and fractures to the outlet during year 35001 and the 38500th year. After that the clean water starts to infiltrate into the whole domain (Figure 34).Conductivity distribution and adaptive triangular mesh for Case 7.
(a)
Conductivity distribution(b)
Adaptive triangular meshPressure distribution, streamlines field, and Average effluent concentration for Case 7.
(a)
Pressure distribution(b)
Streamlines field(c)
Average effluent concentrationConcentration distribution at different time for the first phase for Case 7.
(a)
The 3500th year(b)
The 14000th year(c)
The 35000th yearConcentration distribution at different time for the second phase for Case 7.
(a)
The 38500th year(b)
The 49000th year(c)
The 70000 year


## 4.3. Analysis of the Concentration on the Outflow Boundary at Different Time

We get the effluent concentration on the outlet for every case at different time (Figures8, 12, 16, 20, 24, 28, and 32). For Case 1 to Case 7, they have a common feature that the effluent concentration quickly becomes higher in shorter time than that in the base case during the first phase, and then the concentration increases gradually until it gets to the highest (1.0). This suggests that there are time multiscale phenomena in the contaminant transport in karstic aquifers.  A main reason may be that the vugs and fractures have much higher porosity than that in matrix porous media, and the velocity in the vugs and fractures is also higher than that in matrix porous media. Once the clean water is injected into the inlet, the fractured vuggy porous media, the effluent concentration becomes lower in shorter time than that in the base case except for the Case 4. It is easily to understand that the vugs and fractures are the main pathways for contaminant transport. But it takes longer time for the vuggy fractured porous media to get the lowest concentration than that in base case (Figures 8, 12, 16, 20, 24, 28, and 32). It may be because when the clean water is injected into the inlet, it firstly infiltrates into the vugs and fractures, but there remains much contaminant in the matrix porous media by adsorption. So we can conclude that karst groundwater becomes polluted more easily and in shorter time periods than that in non-karstic aquifers. Even though the aquifers are injected clean water which has been polluted by some contaminants, it will take a much longer time to render the groundwater to be clean again than that in non-karstic areas, because of the significant difference of adsorption between the vugs, fractures, and matrix porous media.


## 5. Conclusions

This paper begins with presenting a mathematical model for some contaminant transport in the fractured vuggy porous media. Two phases are numerically simulated for a process of contaminant and clean water infiltrated in the fractured vuggy porous media by coupling mixed finite element (MFE) method and finite volume method (FVM), both of which are locally conservative, to approximate the model. A hybrid mixed finite element (HMFE) method is applied to approximate the velocity field for the model. The convection and diffusion terms are approached by FVM and the standard MFE, respectively. The pressure distribution and temporary evolution of the concentration profiles are obtained for two phases. The average effluent concentration on the outflow boundary is obtained at different time that shows some obviously different features from the matrix porous media, and the time multi-scale of the effluent concentration on the outlet is observed. The results show how the different distribution of the vugs and the fractures impact the contaminant transport and the effluent concentration on the outlet.



---

*Source: 498098-2011-03-24.xml*