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.
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 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 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 the 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 coefficient , 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 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 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 () dimensional elements, the discrete fracture or vug model in this paper still uses physically sense 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 Here 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, is the gravity acceleration, and are the kinematic (s/m2) and dynamic viscosities of the fluid (kg/(ms)), respectively; and are the density of the fluid (kg/m3) and the (absolute) permeability (m2) of the porous media, respectively. The unknowns are (the pressure head of the fluid mixture (N/m2) and (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 is uniformly symmetric positive definite and bounded. For simplicity of computation, the computational domain of interest is assumed to be polygonal and bounded in with boundary (, stand for the Dirichlet boundary and Neumann boundary, resp.). We take the boundary conditions as follows: where and are the given pressure on and the normal velocity component on ( and can be measured in practice in some points on boundary), respectively, and is the outward norm vector towards.
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: In the system (2.4), and 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. represents a parameter of the distribution coefficient of the considered contaminant species between the matrix and fluid. , the dispersion-diffusion tensor, contributes from molecular diffusion and mechanical dispersion, and it can be computed by where denotes the projected tensor onto the direction, which is calculated by where represents the Euclidean norm of ; (strictly positive) is the molecular diffusivity, and are the transverse dispersivity and the longitudinal dispersivity, respectively, and both are nonnegative. The unknowns as well as 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 is calculated for the matrix and for the vugs, respectively, as follows:
-in the vugs and/or fractures: , , ,
-in the matrix: , .
The boundary conditions for transport system are given by where is the prescribed given concentration on the boundary; denotes the inflow boundary, and denotes the outflow boundary, defined by The initial condition is taken as where is a prescribed given function for the concentration at , 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 inner product over a domain for scalar functions, and denotes inner product for vector functions. Some necessary spaces are given by
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 find and such that where denotes the velocity extension which lets its component agrees with on . 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 th order RT space for the two-dimensional triangular element is given by where represents the direct sum. We employ the case in this paper. Therefore, has the form Restricted to the element, is the polynomial space of degree less than or equal to . We apply space in our coming numerical examples. So our MFE method for approaching the weak formulation for the flow equation is to find and such that 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 space, (3.5) can be expressed as an algebraic linear system with the unknowns (the pressure edge averages) 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 solution and the diffusive flux solution such that where 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 and such that Then we specify a fully discretized algorithm for the transport. We discretize the simulation time into subintervals:. Let ,. Provided that there exists a constant such that , 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 , ?? such that We employ the standard MFE algorithm to solve reactive transport equation (3.7). We here also use 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.
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 concentration and boundary condition in the inflow boundary , . In the first phase, the polluted water infiltration process, with initial concentration in the whole region and boundary condition in the inflow boundary??. In the second phase, again the clean water infiltration process, with initial concentration in the whole region and boundary condition in the inflow boundary , , where is the concentration distribution at in the first phase. The domain considered for all cases in this paper is a bounded rectangular domain of 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 (?m) and the bottom boundaries (?m). The inflow boundary and outflow boundary are taken as the left boundary (?m) and the right boundary (?m) of the domain, respectively. A constant pressure of 0 (the gauge pressure against a reference pressure) is specified on the outflow boundary (?m). In the first phase of the simulation, the polluted water is continuously injected on the left boundary (?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.
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 . 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 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.
(a) Conductivity distribution
(b) Adaptive triangular mesh
(a) Pressure distribution
(b) Streamlines field
(a) The 3500th year
(b) The 14000th year
(c) The 35000th year
(a) The 12000 year
(b) The 15000 year
(c) The 20000 year
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.
(a) Conductivity distribution
(b) Adaptive triangular mesh
(a) Pressure distribution
(b) Streamlines field
(c) Average effluent concentration
(a) The 3500th year
(b) The 14000 year
(c) The 35000th year
(a) The 38500th year
(b) The 49000th year
(c) The 70000th year
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.
(a) Conductivity distribution
(b) Adaptive triangular mesh
(a) Pressure distribution
(b) Streamlines field
(c) Average effluent concentration
(a) The 3500th year
(b) The 14000th year
(c) The 35000th year
(a) The 38500th year
(b) The 49000th year
(c) The 70000th year
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).
(a) Conductivity distribution
(b) Adaptive triangular mesh
(a) Pressure distribution
(b) Streamlines field
(c) Average effluent concentration
(a) The 3500th year
(b) The 14000th year
(c) The 35000th year
(a) The 38500th year
(b) The 49000th year
(c) The 70000th year
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).
(a) Conductivity distribution
(b) Adaptive triangular mesh
(a) Pressure distribution
(b) Streamlines field
(c) Average effluent concentration
(a) The 3500th year
(b) The 14000th year
(c) The 35000th year
(a) The 38500th year
(b) The 49000th year
(c) The 70000th year
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).
(a) Conductivity distribution
(b) Adaptive triangular mesh
(a) Pressure distribution
(b) Streamlines field
(c) Average effluent concentration
(a) The 3500th year
(b) The 14000th year
(c) The 35000th year
(a) The 38500th year
(b) The 49000th year
(c) The 70000 year
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.
(a) Conductivity distribution
(b) Adaptive triangular mesh
(a) Pressure distribution
(b) Streamlines field
(c) Average effluent concentration
(a) The 3500th year
(b) The 14000th year
(c) The 35000th year
(a) The 38500th year
(b) The 49000th year
(c) The 70000th year
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).
(a) Conductivity distribution
(b) Adaptive triangular mesh
(a) Pressure distribution
(b) Streamlines field
(c) Average effluent concentration
(a) The 3500th year
(b) The 14000th year
(c) The 35000th year
(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 (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.
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.
Appendices
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, and are scalar functions or vector functions, which depends on the physical meaning of the quantities and . And the boundary of is denoted by , representing the Dirichlet boundary, representing the Neumann boundary. And the divergence theorem is given by where is the outward unit normal vector towards .
The flow equation considered in this paper here is the same where , , , and is the velocity extension that make the normal component of agreeing with on , where is the prescribed given normal component on .
We begin with rewriting the first equation in (A.2) as the following one: Multiplying (A.3) by and integrating over , we have For (A.5), by the divergence theorem and the definition of , we can obtain By the definition of inner product, we can get where is the value of on .
We rewrite (A.7) with the notation of inner product as follows: Because , we can get Multiplying the second one in the flow equation by , similarly, we can obtain Combining (A.9) and (A.10), we get the weak formulation of the flow equation
B. Matrix Formulation of Flow Equation
We derive the matrix formulation of the flow equation. We begin with a single triangular element . By space, in (3.5) can be written as where is a shape function of the space and represents the total flux across an edge of element . Taking as the test function, the total flux is given by By (B.2), we have where , , , , . As a result, the flow in (B.2) passing through each edge is given by a function of the cell pressure average and edge pressure average . For simplicity, we rewrite the equation as where , .
For the second one in (3.5), combining (B.1), we can obtain where , ().
Substituting (B.3) into (B.5) we can get Hence by the continuity of the flux across the interelement boundaries, then is given by Thus by (B.2), (B.7), and (B.8), we can obtain the algebraic linear system of (the pressure average) as follows: where , , , , , . and stand for the number of elements and the number of edges not belonging to in the mesh, respectively. represents a vector of size 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 by MFE method in (2.7), and we rewrite (C.1) to get Now we can get Multiplying the first equation in (C.3) by , the second by , and the last by , we can obtain Using the upwind value of the concentration, the boundary condition, and the initial condition, we now get the weak formulation
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 space, the in (3.8) can written as where is a basis function of space, is the total diffusive flux across the edge . We take the test function as and integrate over , using the Green formulation, and therefore the total diffusive flux is given by where , , and . So the diffusion term in the first equation in (3.8) and we can get where .
By (D.2) and (D.4), we have the ordinary differential equations of the concentration cell averages and the diffusive flux as follows: For the advection component in system (3.8), Because is a function of space, . The given equation above is written as Denote the velocity across each edge by , and then give where is the outward unit normal vector to each edge. Denote the edge-element signed adjacency matrices by the entry of which equals to 1 if for each edge, and the entry of equals to -1 if . And therefore the upwind concentration is given by or . Let , , and , using the boundary condition, and then we can obtain the flux across the edges as follows: where and represent the edge-boundary adjacency matrix and the concentration on the inflow boundary, respectively. We denote edge-adjacency matrix by . And then the total divergence amount is given by Equation (D.10) gives an ordinary differential equations as follows: Now we finally can obtain the ordinary differential equations by (D.4) and (D.11) for the transport system where , , , .
Acknowledgments
This paper is supported by the Foundation for Young Talents of Guizhou Province under Grant Z073245. The authors thank its support. The authors also take this opportunity to thank the Computational Transport Phenomena Laboratory (CPTL), Division of Physical Sciences and Engineering (PSE), King Abdullah University of Science and Technology (KAUST) for supporting this paper.