Numerical Simulation of Pollutant Transport in Fractured Vuggy Porous Karstic Aquifers

This paper begins with presenting a mathematical model for contaminant transport in the fractured vuggy porous media of a species of contaminant (cid:2) PCP (cid:3) . Two phases are numerically simulated for a process of contaminant and clean water inﬁltrated in the fractured vuggy porous media by coupling mixed ﬁnite element (cid:2) MFE (cid:3) method and ﬁnite volume method (cid:2) FVM (cid:3) , both of which are locally conservative, to approximate the model. A hybrid mixed ﬁnite element (cid:2) HMFE (cid:3) method is applied to approximate the velocity ﬁeld for the model. The convection and di ﬀ usion terms are approached by FVM and the standard MFE, respectively. The pressure distribution and temporary evolution of the concentration proﬁles are obtained for two phases. The average e ﬄ uent concentration on the outﬂow boundary is obtained at di ﬀ erent time and shows some di ﬀ erent features from the matrix porous media. The temporal multiscale phenomena of the e ﬄ uent concentration on the outlet are observed. The results show how the di ﬀ erent distribution of the vugs and the fractures impacts on the contaminant transport and the e ﬄ uent concentration on the outlet. This paper sheds light on certain features of karstic groundwater are obtained.


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 Figure 1.
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 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 the x − 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 Darcyweighted 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 coefficient K d , 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 T 1 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.

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.

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 u K∇p 0,

2.1
Here K denotes the hydraulic conductivity i.e., the proportionality constant for the flow of water through a porous media, m/s as follows: In the two equations above, g is the gravity acceleration, ν and η are the kinematic s/m 2 and dynamic viscosities of the fluid kg/ m•s , respectively; ρ and κ are the density of the fluid kg/m 3 and the absolute permeability m 2 of the porous media, respectively.The unknowns are p the pressure head of the fluid mixture N/m 2 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 R 2 with boundary Γ ∂Ω Γ D ∪ Γ N Γ D , Γ N stand for the Dirichlet boundary and Neumann boundary, resp. .We take the boundary conditions as follows: where p B and u B are the given pressure on Γ D and the normal velocity component on Γ N p B and u B can be measured in practice in some points on boundary , respectively, and n is the outward norm vector towards u.

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 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.K d 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 where E u denotes the projected tensor onto the u direction, which is calculated by where |u| u 2 1 u 2 2 1/2 represents the Euclidean norm of u u 1 , u 2 ∈ R 2 ; d m strictly positive is the molecular diffusivity, d t and d l are the transverse dispersivity and the longitudinal dispersivity, respectively, and both are nonnegative.The unknowns c f as well as c s 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 where φ e φ ρK d is calculated for the matrix and for the vugs, respectively, as follows: -in the vugs and/or fractures: φ 1, K d 0, φ e φ 1.0, -in the matrix: φ φ m , φ e φ ρK d .The boundary conditions for transport system are given by where c B is the prescribed given concentration on the boundary; Γ inflow denotes the inflow boundary, and Γ outflow denotes the outflow boundary, defined by The initial condition is taken as where c 0 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.

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.

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.

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 L 2 Ω inner product over a domain Ω ∈ R 2 for scalar functions, and L 2 Ω 2 denotes inner product for vector functions.
Some necessary spaces are given by

3.1
Based on these spaces, we give the weak formulations for flow equation and reactive equation as follows, respectively.

Weak Formulation
For the flow equation, the weak formulation is to find p ∈ V and u ∈ W 0 N E u B such that where E u B denotes the velocity extension which lets its component agrees with u B on Γ N .There is the derivation of 3.2 in Appendix A.

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 where ⊕ represents the direct sum.We employ the case r 0 in this paper.Therefore, W h K has the form Restricted to the element, P r K is the polynomial space of degree less than or equal to r.We apply RT 0 space in our coming numerical examples.So our MFE method for approaching the weak formulation for the flow equation is to find

3.5
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 the RT 0 space, 3.5 can be expressed as an algebraic linear system with the unknowns P E the pressure edge averages The derivation of 3.6 is presented in Appendix B.

Weak Formulation
The weak formulation for the reactive transport system is to find the concentration solution c ∈ V and the diffusive flux solution w ∈ W 0 such that where c * f represents the upwind value of the concentration on an edge.There is the derivation of 3.7 in Appendix C.

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 c fh ∈ V h and w h ∈ W 0 h such that Then we specify a fully discretized algorithm for the transport.We discretize the simulation time 0, T into n subintervals: 0 Provided that there exists a constant C > 0 such that Δt ≤ C min l Δt l , 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 find

3.9
We employ the standard MFE algorithm to solve reactive transport equation 3.7 .We here also use RT 0 on the reactive transport system, and thus the linear system 2.4 can be specified as the following algebraic linear ordinary differential equation system Ordinary differential equations 3.10 are solved by backward Euler method in this paper.The derivation of 3.10 is presented in Appendix D.

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 concentration c 1 x, 0 0 and boundary condition in the inflow boundary c 1 x, t 0, x ∈ Γ inflow .In the first phase, the polluted water infiltration process, with initial concentration in the whole region c 2 x, 0 0.0 and boundary condition in the inflow boundaryc 1 x, t 1.0, x ∈ Γ inflow .In the second phase, again the clean water infiltration process, with initial concentration in the whole region c 3 x, 0 c 1 x, T 1 and boundary condition in the inflow boundary c 3 x, t 0, x ∈ Γ inflow , where c 1 x, T 1 is the concentration distribution at T 1 in the first phase.The domain considered for all cases in this paper is a bounded rectangular domain Ω ⊂ R 2 of 0, 0.80 m × 0, 0.80 m 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 .
Case 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.

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, T 1 .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 RT 0 − 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.

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

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

Case 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 Figure 15.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 .

Case 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 Figure 19.
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 .

Case 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 Figure 23.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 .

Case 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 Figure 27.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.

Case 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 Figure 31.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 .

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 Figures 8,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.

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.

A. Weak Formulation of Flow Equation
We represent the derivation of the weak formulation of the flow equation here.
We at first denote the inner product notation by the following formulation where Ω is a plane area, u and v are scalar functions or vector functions, which depends on the physical meaning of the quantities u and v.And the boundary of Ω is denoted by      where n is the outward unit normal vector towards Γ.
The flow equation considered in this paper here is the same Ω and E u B is the velocity extension that make the normal component of u agreeing with u B on Γ N , where u B is the prescribed given normal component on Γ N .We begin with rewriting the first equation in A.2 as the following one:      By the definition of inner product, we can get where p B is the value of p on Γ.We rewrite A.7 with the notation of inner product as follows: Because w ∈ W 0 N , we can get A.9  Multiplying the second one in the flow equation by v ∈ V , similarly, we can obtain Combining A.9 and A.10 , we get the weak formulation of the flow equation A.11  By B.2 , we have As a result, the flow u K,E in B.2 passing through each edge is given by a function of the cell pressure average p K and edge pressure average p K,E .For simplicity, we rewrite the equation as where  by the continuity of the flux across the interelement boundaries, then u K,E is given by

B.8
Thus by B.2 , B.7 , and B.8 , we can obtain the algebraic linear system of P E the pressure average as follows: A T P E J, B.9 where , E / ∈ Γ D .N K and N E stand for the number of elements and the number of edges not belonging to Γ D in the mesh, respectively.J represents a vector of size N E denoting the boundary conditions.

C. Weak Formulation of the Transport Equation
We derive the weak formulation of the transport equation in this appendix.We let

D. MFE Method to Solve the Transport System
We apply the standard MFE method to solve the transport system 3.8 .We handle the second one in the system similarly to handling flow equation.Based on RT 0 space, the w h in 3.8 can written as where w K,E is a basis function of RT 0 space, w K,E is the total diffusive flux across the edge E. We take the test function w as w K,E and integrate over Ω, using the Green formulation, and therefore the total diffusive flux is given by  where n i is the outward unit normal vector to each edge.Denote the edge-element signed adjacency matrices by F the entry of which equals to 1 if u ⊥ iE > 0 for each edge, and the entry of F − equals to −1 if u ⊥ iE < 0. And therefore the upwind concentration is given by F C or F − C. Let u ⊥ u ⊥ iE N E , u ⊥ max 0, u ⊥ , and u ⊥ − min 0, u ⊥ , using the boundary condition, and then we can obtain the flux across the edges as follows: where B and C B represent the edge-boundary adjacency matrix and the concentration on the inflow boundary, respectively.We denote edge-adjacency matrix by F F F − .And then the total divergence amount is given by Equation D.10 gives an ordinary differential equations as follows:

Figure 2 :
Figure 2: Injection on the inflow boundary for the two phases of the simulation.

Figure 3 :
Figure 3: Conductivity distribution and adaptive triangular mesh for Case 0.

Figure 4 :
Figure 4: Pressure distribution and streamlines field for base case.

Figure 5 :
Figure 5: Concentration distribution at different time for the first phase for base case.

Figure 6 :
Figure 6: Concentration distribution at different time for the second phase for base case.

Figure 7 :
Figure 7: Conductivity distribution and adaptive triangular mesh for Case 1.

Figure 8 :
Figure 8: Pressure distribution, streamlines field, and average effluent concentration for Case 1.

Figure 9 :
Figure 9: Concentration distribution at different time for the first phase for Case 1.

Figure 10 :
Figure 10: Concentration distribution at different time for the second phase for Case 1.

b
Adaptive triangular mesh

Figure 11 :
Figure 11: Conductivity distribution and adaptive triangular mesh for Case 2.

Figure 12 :
Figure 12: Pressure distribution, streamlines field, and average effluent concentration for Case 2.

Figure 13 :
Figure 13: Concentration distribution at different time for the first phase for Case 2.

Figure 14 :
Figure 14: Concentration distribution at different time for the second phase for Case 2.

b
Adaptive triangular mesh

Figure 15 :
Figure 15: Conductivity distribution and adaptive triangular mesh for Case 3.

Figure 16 :
Figure 16: Pressure distribution, streamlines field, and average effluent concentration for Case 3.

Figure 17 :
Figure 17: Concentration distribution at different time for the first phase for Case 3.

Figure 18 :
Figure 18: Concentration distribution at different time for the second phase for Case 3.

b
Adaptive triangular mesh

Figure 19 :
Figure 19: Conductivity distribution and adaptive triangular mesh for Case 4.

Figure 20 :
Figure 20: Pressure distribution, streamlines field, and average effluent concentration for Case 4.

Figure 21 :
Figure 21: Concentration distribution at different time for the first phase for Case 4.

Figure 22 :
Figure 22: Concentration distribution at different time for the second phase for Case 4.

b
Adaptive triangular mesh

Figure 23 :
Figure 23: Conductivity distribution and adaptive triangular mesh for Case 5.

Figure 24 :
Figure 24: Pressure distribution, streamlines field, and average effluent concentration for Case 5.

Figure 25 :
Figure 25: Concentration distribution at different time for the first phase for Case 5.

Figure 26 :
Figure 26: Concentration distribution at different time for the second phase for Case 5.

b
Adaptive triangular mesh

Figure 27 :
Figure 27: Conductivity distribution and adaptive triangular mesh for Case 6.

Figure 28 :
Figure 28: Pressure distribution, streamlines field, and average effluent concentration for Case 6.

Figure 29 :
Figure 29: Concentration distribution at different time for the first phase for Case 6.

Figure 33 :
Figure 33: Concentration distribution at different time for the first phase for Case 7.

Figure 34 :
Figure 34: Concentration distribution at different time for the second phase for Case 7.

w
−D u ∇c, C.1 by MFE method in 2.7 , and we rewrite C.1 to get∇c −D −1 u w. n c B u • n, x ∈ Γ inflow , −w • n 0, x ∈ Γ outflow , c f x, 0 c 0 x , x ∈ Ω. C.3Multiplying the first equation in C.3 by w, the second by v, and the last by w, we can obtain∂ φ e c f ∂t , v ∇ • uc f w , v 0, v , v ∈ V, t ∈ 0, T , ∇c, w −D −1 u w, w , w ∈ W 0 , c f x, 0 , v c 0 x , v , v ∈ V, t 0, uc w • n c B u • n, x ∈ Γ inflow , −w • n 0, x ∈ Γ outflow .C.4Using the upwind value of the concentration, the boundary condition, and the initial condition, we now get the weak formulation ∈ V, t ∈ 0, T , −D −1 w, w c f , ∇ • w 0, ∀w ∈ W 0 , t ∈ 0, T , c f , v c f,0 , v , ∀v ∈ V, t 0. C.5 v K • v K K φ e dK v K 1 .By D.2 and D.4 , we have the ordinary differential equations of the concentration cell averages C and the diffusive flux W K as follows: Because v is a function of RT 0 space, ∇v 0. The given equation above is written as across each edge by u i , and then give u ⊥ iE u i • n i , i 1, 2, . . ., N E , D.8

Table 1 :
Standard parameters for the mathematical model in practice.