Numerical Simulation of Groundwater Flow and Solute Transport in a Karst Aquifer with Conduits

Abstract

Numerical simulation of groundwater flow and solute transport in a dual-permeability karst aquifer is a challenging issue, since groundwater flow in karst conduit network can become non-darcian even turbulent flow. The discrete-continuum model is a relatively new modeling method, which accounts for turbulent and laminar flow in karst aquifer. MODFLOW-CFP (Conduit Flow Process) is compared to the MODFLOW, a numerical code based on Darcy law, to evaluate the accuracy in a sub-regional scale karst aquifer. MODFLOW-CFP is more accurate than the MODFLOW when comparing the head simulation results with field measurements. After that, the CFPv2 and UMT3D numerical models are applied in the WKP to establish a sub-regional scale model to simulate chloride transport processes in the last four decades, and to predict contamination development. Numerical simulation results indicate sprayfields are the major chloride source in the study region. Conduit networks significantly control solute transport and contaminant distribution in the study region. Chloride transports through conduits rapidly and spread to several large contamination plumes in a short period. Chloride concentration started to increase in 1980s due to the operation of sparyfield. Solute transport simulation results by discrete-continuum models are more accurate because of the precise description of conduit network.

1. Introduction

Karst aquifer systems underlie approximately 10–20% of the Earth's landmass and supply potable water to nearly 25% of the world's population [1]. In the United States, karst carbonate aquifers supply almost 52% of all bedrock aquifer withdrawals on a yearly basis [2, 3]. Karst aquifers are typically characterized by relatively large void spaces and loose porous media, which make the karst aquifers to be the most productive aquifer systems in the world, for example, the Floridan aquifer in the USA [4]. The open and porous nature of karst aquifers, combined with the dissolution of joints and fractures within carbonates over time, generate complex subsurface conduit systems and dual‐permeability flow regimes [5]. Groundwater flow and solute transport processes in conduit networks are generally more rapid than that in the surrounding carbonate porous media due to significantly higher hydraulic conductivity and porosity in the conduit system [6, 7]. As a result, the water and solute interchanges between the conduit and matrix are particularly important in a dual‐permeability karst aquifer system. Karst aquifers are particularly vulnerable when one considers the problems of groundwater contamination [8]. During high‐flow events, contaminants in the conduit system can be rapidly transported through the aquifer or actively pushed into the carbonate matrix when the water pressure is high in the conduit system. The process is reversed during low‐flow events, and contaminants are slowly released from the surrounding matrix into the conduit system [9]. Therefore, groundwater contamination and seawater intrusion in a karst aquifer can persist for a long time because of the retention and release effect between the conduits and matrix domains [10, 11].

Laboratory experiments using physical models with artificial gypsum or sandbox analog were taken to simulate karst groundwater flow and solute transport [12]. Li [13] evaluated the solute transport and retention in a karst aquifer using two categories of laboratory experiments. Faulkner et al. [14] designed a sandbox laboratory experiment to simulate the interaction between conduit and porous medium domain. In addition to laboratory studies, many numerical models have been developed to study groundwater flow and solute transport in dual‐permeability systems, as well as chemical reaction and carbonate dissolution in karst aquifers [15]. For example, a limestone dissolution continuum model coupled with conduit flow was proposed to simulate groundwater flow and karst evolution [16]. Lauritzen et al. [17] coupled laminar flow and carbonate dissolution in two‐dimensional (2D) pipe networks to study groundwater flow and karst development in a limestone aquifer. Groves and Howard [18] used 2D pipe networks to simulate conduit development processes under laminar flow conditions at field scales, and the simulation method was later extended to turbulent flow by Howard and Groves [19]. Kaufmann and Braun [20] coupled a pipe network with a continuum system to study karst development and indicated that early karstification might be enhanced by the presence of a diffuse flow system.

The coupling of nonlinear or turbulent flow within the conduit network and Darcian laminar flow in the porous medium is a challenging issue in numerical simulation of a dual‐permeability karst aquifer [5, 21, 22]. Darcy equation of continuum groundwater flow is applicable for laminar flow in the porous medium that is linear to hydraulic gradient but is not accurate for nonlaminar or turbulent flow in the conduit. Some hybrid discrete‐continuum numerical models were developed to couple Darcian flow in the porous medium with nonlaminar channel flow in the conduit, such as carbonate aquifer void evolution (CAVE) [23, 24], MODFLOW‐CFPM1 [21] and CFPv2 [25, 26]. Those discrete‐continuum models were also applied and evaluated in a number of studies [22, 27–30].

The solute transport simulation in the discrete‐continuum model involves the coupling of 1D solute transport within the conduit network and the 2D/3D solute transport in the porous medium domain. MT3DMS is a widely used modular 3D model to simulate solute transport in a porous medium, which is based on the solution of groundwater flow from MODFLOW simulation [31]. Therefore, MT3DMS is only applicable for solute transport in porous medium aquifer [31]. Spiessl [32] and Spiessl et al. [33] modified the source code of MT3DMS and developed the reactive hybrid transport RUMT3D model. RUMT3D was coupled with CAVE to calculate groundwater flow in the conduits and then solve solute transport equations in the conduit and porous medium domains of a karst aquifer. Reimann et al. [34] enabled MODFLOW‐CFP simulation of unsaturated flow in conduits and water exchange between matrix and variably filled conduits. Furthermore, conduit‐associated drainable storage (CADS) was added in MODFLOW‐CFPM1 and fixed head‐limited flow boundary condition in Reimann et al. [25]. Based on those studies, a research version of CFPv2 has been recently developed by Reimann et al. [26] and used in Xu et al. [29] to simulate the flow of laminar and nonlaminar, as well as nitrate‐N transport in conduits and matrix with an updated UMT3D. The CFPv2 and UMT3D models were combined to calculate solute advective exchange between conduit and medium, as well as transport process within a porous medium by advection and dispersion.

2. Study site and model setup

The Woodville Karst Plain (WKP) is a geographic area in the south of Tallahassee, Florida, characterized by prominent karst features such as sinking streams, sink holes, submerged caves systems and springs (Figure 1). The upper Floridan aquifer is the primary aquifer, which is composed of the Oldsmar to Chattahoochee/St. Marks formations, with the Clayton formation on the bottom and Hawthorn group acting as the confining unit on the top. Groundwater flow is generally from north to south in the Floridan aquifer. The impermeable Hawthorn group is thin or even absent in the Woodville Karst Plain. Rainwater is able to infiltrate easily and quickly, dissolve the limestone structure and then form substantial networks of karst conduits and sinkholes with a substantial amount of water flows. There are enormous springs existed within the WKP, and the three primary spring discharge points for the cave networks are Wakulla Spring, Spring Creek Springs Group and St Marks Spring. Figure 1 shows the study area and the hydraulic head contours based on field measurements.

A subregional scale MODFLOW‐2000 and MT3DMS models of the Woodville Karst Plain is created by Davis et al. [35], which simulated long‐term groundwater flow and nitrate transport in the Woodville Karst Plain. The MODFLOW‐2000 model by Davis et al. [35] was created based on a previous regional‐scale model, which includes the entire north‐central Florida and southeast Georgia by Davis and Katz [36]. In the study by Davis et al. [35], large hydraulic conductivity cells are used in the cells of “cave network” to account for flow in submerged conduits within the Woodville Karst Plain. The simulated heads, flow and transport within the Woodville Karst Plain by Davis et al. [35] match with the observations very well.

In Gallegos et al. [22], the CFPM1 model of the Woodville Karst Plain is based on the MODFLOW‐2000 model by Davis et al. [35] and thus very similar from a conceptual standpoint (Figure 2). The conceptual model is divided into two layers, which represent the upper Floridan aquifer. Below the two layers are low permeability sediments, which act as a no‐flow boundary. The north boundary and south boundary are set as specified head and constant head boundary, respectively. On the other hand, the southeastern boundary acts as a no‐flow boundary. The formation of sinkholes, conduits and caves in the Floridan aquifer is due to the dissolution of limestone with water flowing through these karst features. Sand and clay are observed to cover the upper portion of the aquifer (layer 1) as a low hydraulic conductivity unit. The lower portion of the aquifer (layer 2) has relatively higher hydraulic conductivities with little to no infilling. Therefore, most of the karst conduit networks are located in layer 2. In this study, the Floridan aquifer is simulated as a confined aquifer. Discharge is only simulated at outlets of conduit networks as springs, and at the south boundary that represents the shoreline of Gulf of Mexico. Only the discharges of three springs are simulated in the numerical model: Wakulla Spring, Spring Creek Springs Group and St Marks Spring.

There are 288 rows and 258 columns in the discretization of grid of the subregional CFPM1 model. Model cells are 500 by 500 ft horizontally. The top and bottom of numerical mode are determined from the regional‐scale model by Davis and Katz [36]. There are two layers in the CFPM1 model, according to the conceptual model. The upper portion of the upper Floridan aquifer in layer 1 has relatively low hydraulic conductivities. The lower portion of the upper Floridan aquifer in layer 2 has relatively large hydraulic conductivities. The lateral boundaries were set as specified head boundaries in both two layers, except the southeastern portion of the model is set as no‐flow boundary. The bottom of lower layer is also considered as a no‐flow boundary. Drain package is used to simulate the rivers [37]. Hydraulic conductivity values for the matrix were originally based on the calibrated values from Davis et al. [35]. There are five sources of recharge in the model, including the inflow across lateral boundaries of the model, net precipitation, creek flow into sinkholes, irrigation at sprayfields and discharges from onsite sewage disposal systems (e.g., septic tanks).

In Davis et al. [35], grid cells of submerged conduits and caves are simulated as large values of hydraulic conductivity. The actual mapped caves were only small portions of the high conductivity cells in the model. The remaining very high conductivity cells were placed in their respective locations based on inference from field observations that indicated probable cave locations, which are indicated by tracer tests and the presence of numerous sinkholes due to heavy dissolution in the subsurface. The high tannic flows at lost Creek Sink resulting in flows of dark water at Spring Creek Springs Group also indicate indicating a hydraulic connection between Wakulla Springs and the Spring Creek Springs Group [35, 38, 39]. Cave diver observations found that the water was flowing southward at the southernmost tip of the explored Wakulla caves system.

Figure 1.

Subregional model boundary of the Woodville Karst Plain. Modified from Davis et al. [35].

Figure 2.

Conceptual model of the Woodville Karst Plain subregional model. Modified from Davis et al. [35].

Figure 3.

Pipe network in the CFPM1 subregional model of the Woodville Karst Plain. Modified from Gallegos et al. [22].

There are 37 miles explored and mapped conduit networks by cave divers that connect to Wakulla spring [35, 38], which is shown as the brown lines in Figure 1. However, the conduit network distribution in the WKP is extended and connected with numerous springs and sinks, including Wakulla Spring, Spring Creeks Springs, Lost Sink, etc. In Gallegos et al. [22], a pipe network showing in Figure 3 was created using CFPM1 to replace the very high hydraulic conductivity cells from the study of Davis et al. [35]. For the sake of comparison between the CFPM1 and MODFLOW‐2000 models, the discrete pipes and nodes of conduit were placed in the exact same cells in layer 2 as the cells with large hydraulic conductivity values in the original MODFLOW‐2000 model. In CFPM1 model, conductivity values of the high hydraulic conductivity cells were assigned to be the same as the surrounding hydraulic conductivity cells. There are 1083 pipe nodes and 1087 pipe segments of the conduit network system in the numerical model. In the discrete model, the conduits were subdivided into discrete pipe network groups, based on the locations of recharge or discharge points such as sinks or springs. Discrete conduit network were assumed to be located at the horizontal and vertical center of the model cells; therefore, pipe nodes were assigned to the center of the cells. Cave divers measured the actual diameters of conduit network within the Wakulla Cave System; however, all other diameters for pipes were estimated during calibration or estimated based on the available field data. Tortuosity was calculated for each major cave segment mapped by the cave divers. The mean height of the pipe wall microtopography (i.e., wall roughness) was assumed to be 0.10. Three major springs, including Wakulla Spring, Spring Creek Springs Group and St. Marks Spring, are set as three pipe nodes of constant head in layer 1.

In the study of Xu et al. [29], effective porosity governs solute advection and dispersion transport simulation. In groundwater flow simulation, MODFLOW calculates groundwater velocity within model cell and assumes that groundwater fills the entire cells. However, groundwater velocity is calculated at a more realistic rate by introducing the effective porosity, which restricts the groundwater movement to the percentage of the cell. Effective porosity is difficult to assess accurately, especially in karst terrains. The simulated effective porosity for layer 1 is set as uniformly 0.01. An effective porosity of 0.003 in the area was used in entire model cells in layer 2. However, effective porosity is set as 0.03 in area close to the SEF sprayfield, because aquifers there are quite loose and small conduits are so complicated that we cannot only assume the conduit nodes and tubes using several big tubes by discrete pipes to calculate groundwater flow and solute transport.

3. Numerical simulation of groundwater flow in a karst aquifer

3.1. Governing equations of groundwater flow

According to Darcy equation and mass conservation, the 3D continuum equation for groundwater flow in a porous medium including source/sink term is described as follows:

$\frac{\partial }{\partial x}\left(K_{xx}\frac{\partial h}{\partial x}\right)+\frac{\partial }{\partial y}\left(K_{yy}\frac{\partial h}{\partial y}\right)+\frac{\partial }{\partial z}\left(K_{zz}\frac{\partial h}{\partial z}\right)\pm W= S_s\left(\frac{\partial h}{\partial t}\right)$E1

where $K_{\text{xx}},{\text{K}}_{\text{yy}}$ and ${\text{K}}_{\text{z}\text{z}}$ are values for hydraulic conductivity [LT^-1] along the x, y, z axis, respectively; $h$ is the hydraulic head [L]; $W$ is the volumetric flux per unit volume [T^-1] as sink/source term; $S_s$ is the specific storage [L^-1].

Darcy equation is only appropriate for laminar flow with Reynolds number less than 10 but not for nonlaminar or turbulent flow in the conduit [40]. Therefore, Darcy‐Weisbach equation is used to simulate groundwater flow in the high permeable karst conduit networks in which a conduit is conceptualized as a pipe,

$\Delta h= h_L=f\frac{\Delta l}{d}\frac{V^2}{2g}$E2

where $\Delta h$ or $h_L$ is the head loss [L] measured along the pipe length $\Delta l$ [L]; $f$ is the friction factor [dimensionless]; $d$ is the pipe diameter [L]; $V$ is the mean velocity [LT^-1]; $g$ is the gravitational acceleration constant [LT^-2].

Darcy‐Weisbach equation can be reformulated to solve for volumetric flow rate in the conduit networks by MODFLOW‐CFP [21] and CFPv2 [26]. MODFLOW‐CFP and CFPv2 are able to solve nonsaturated pipe flow but not applied in this study, because all layers are assumed as confined aquifers and conduits are fully saturated. Caves and conduits are defined as a discrete pipe network consisting of cylindrical tubes and nodes within the cells of grid defined in the model. The advective exchanges between porous media and conduit nodes are assumed to be linear with head difference as follows:

$Q_{ex}={\alpha }_{j, i, k}\left(h_n-h_{j,i,k}\right)$E3

where $Q_{ex}$ is the volumetric flow exchange rate [L³T^-1]; ${\alpha }_{j, i, k}$ is the pipe conductance at MODFLOW cell $j,i,k$ [L²T^-1]; $h_n$ is the head [L] at pipe node n located at the center of the MODFLOW cell; and $h_{j,i,k}$ is the head [L] in the encompassing MODFLOW cell $j,i,k$.

3.2. Flow modeling application

The discharges into sinkholes were turned off in the simulation of 1991. The pipe nodes of Wakulla Spring, Spring Creek Springs Group and St Marks Spring are set as constant5, 0 and 9 ft, respectively. The yearly averaged head values are used for the three springs, respectively. In the 2006 simulation, all discharges into sinkholes were turned on. The pipe node for Wakulla Spring, Spring Creek Springs Group and St Marks Spring were set as constant 5, 9.5 and 9 ft, respectively. After 2006, the flow condition of Spring Creek changed, where freshwater discharges stopped and water began flowing northward from Spring Creek to Wakulla Spring. Therefore, the head at Spring Creek Springs Group was increased to 9.5 ft as calculated equivalent freshwater head.

The simulated and observed water levels for both the 1991 and 2006 hydraulic head data sets are presented in Figure 4. With the exception of two points, simulated hydraulic heads were within the calibration criterion of ±5 ft. One point of the 1991 data falls outside the calibration criterion limits. This is caused by problematic well that is located in an area where the hydraulic gradient is very steep. The simulated heads also did not fit well within the calibration criterion in Davis et al. [35]. The point of 2006 data falls outside the calibration criterion is located in the sprayfield in the northeast corner of the model, which may have large uncertainty because of the regional flow from north boundary. The residual mean was 1.03 ft with a residual standard deviation of 3.84 ft for the 1991 hydraulic head data set. On the other hand, the residual mean was 0.09 ft with a residual standard deviation of 2.32 ft for the 2006 data set. Generally speaking, simulation results are able to match reasonably well with the observed heads in monitoring wells.

Figure 4.

Observed versus simulated heads for the calibrated subregional flow model.

4. Numerical simulation of solute transport in karst aquifer

4.1. Governing equations of solute transport

Solute transport modeling in a porous medium needs to solve groundwater flow advection, molecular diffusion and mechanical dispersion. The second‐order partial differential transport equation that applied in MT3DMS [31, 41] and UMT3D [33] has been simplified as 2D transport equation in this study as follows:

$\frac{\partial \left(\theta C\right)}{\partial t}=\frac{\partial }{\partial x_i}\left(\theta D_{ij}\frac{\partial C}{\partial x_j}\right)-\frac{\partial }{\partial x_i}\left(\theta v_{i}C\right)+q_s{C}_s$E4

where $\theta$ is the porosity of the porous medium [dimensionless]; $v_i$ is the seepage or linear pore water velocity [LT^-1], which is related to the specific discharge or Darcian flux through the relationship, $v_i=q_i/\theta$;$C$ is the solute concentration [ML^-3]; $D_{ij}$ is the hydrodynamic dispersion coefficient tensor [L²T^-1];$C_s$ is the solute concentration of water entering from sources or flowing out from sinks [ML^-3];$q_s$ is the volumetric flow rate per unit volume of aquifer representing fluid source (positive) and sink (negative) [T^-1].

Solute transport in the conduit is describe by the 1D advection equation as well, while mechanic dispersion within the conduit is ignored in the VDFST‐CFP model,

$\frac{\partial C_l}{\partial t}=-q_l\frac{\partial C_l}{\partial x}$E5

where $C_l$ is the solute concentration in conduit tube l [ML^-3]; $q_l$ is the conduit flow velocity in conduit tube l [LT^-1], which could be calculated by volumetric flow rate $Q$ from the Darcy‐Weisbach equation [L³T^-1]. It is noted that there are no sink/source terms in the conduit transport equation, because the sink/source terms and exchanges with the surrounding porous media are computed at the conduit nodes only.

Advective exchange rates between a conduit node and the surrounding porous medium cells are determined by the following:

$K_{ex}=\{\begin{array}{c}\frac{Q_{n,ex}^{+}{C}_{i,k}}{V_{i,k}},Q_{n,ex}^{+}>0 \\ \frac{Q_{n,ex}^{+}{C}_n}{V_{i,k}},Q_{n, ex}^{+}<0\end{array}$E6

where $K_{ex,n}$ is the advective exchange rate between a conduit node n and respective matrix cell i, k [ML^-3 T^-1]; $Q_{n,ex}^{+}$ is the exchange flow rate [L³T^-1] at conduit node n ($Q_{n,ex}^{+}>0$, flow direction is from matrix to conduit node; $Q_{n,ex}^{+}<0$, flow direction is from conduit node to matrix); $C_{i,k}$ is the solute concentration of respective cell i, k in the porous medium at conduit node n [ML^-3]; $C_n$ is the nodal concentration at conduit node n [ML^-3]; $V_{i,k}$ is the volume of respective cell i, k of the porous medium at conduit node n [L³].

Spiessl et al. [33] pointed out that mass transport within a conduit networks is determined by the flow velocity within the conduits, exchange coefficients between the conduit nodes and porous matrix, the magnitude of conduit sink/source terms and the lengths of conduit tubes. Mathematically, a weighted arithmetic mean of the nodal concentration value $C_n$ [ML^-3] at conduit node n can be expressed as follows [33]:

$C_n=\frac{{{\displaystyle \sum }}_f^2{Q}_{n,l}^{f+}{C}_{n,l}^f+Q_{n,ex}^{+}{C}_{i,k}+{{\displaystyle \sum }}_s{Q}_{n,s}^{+}{C}_{n,s}}{{{\displaystyle \sum }}_f^2{Q}_{n,l}^{f+}+Q_{n,ex}^{+}+{{\displaystyle \sum }}_s{Q}_{n,s}^{+}}$E7

where superscript f indicates either forward or backward direction of the pipe connected to node n;$C_{n, l}^f$ is the concentration of tube l connected to face f of the conduit node n [ML^-3]; $C_{n, s}$ is the concentration of the source or sink term to the conduit node n [ML^-3].The superscript + represents the inflow terms at conduit node n, which means that only inflow terms are used to compute the nodal concentration; $Q_{n, l }^{f+}$ is the discharge of tube l connected to face f into the respective conduit node n [L³T^-1]; $Q_{n, s}^{+}$ is the volumetric flow rate of a source term at conduit node n [L³T^-1].

4.2. Chloride transport modeling at the WKP

A long‐term (1968–2018) chloride transport process is simulated using the coupled CFPv2 and UMT3D models, based on the MODFLOW‐2000 numerical simulation in the previous studies. Both the simulated chloride concentrations in some monitoring wells and the extents of chloride plume in the WKP are presented as follows.

Chloride level is significantly affected by the application of sprayfields near the City of Tallahassee. Simulation results and water quality measurements at two monitoring wells near the sprayfields are compared. Chloride concentrations also began to increase from background level of less than 5 mg/L after several months of the west part of SEF sprayfield became operational in 1980 for SE‐22 and in 1986 for SE‐21 when the east part of SEF sprayfield became operational. The most severe chloride contaminations at wells SE‐22 (Figure 5, right) can be higher than 35 mg/L in 2000, and chloride concentrations at wells SE‐21 (Figure 5, left) were also as high as 20 mg/L. Chloride levels continued to increase ever since. Chloride concentrations by simulation in those monitoring wells seem to be stabilized at high concentrations after 2005 to now, and also estimate to keep at this value till 2018. Input from the west boundary is a major chloride source, which also indicates very important groundwater inflow from the west.

Figure 5.

Chloride simulation results at SE‐21 and SE‐22 monitoring wells near the Tallahassee SEF sprayfield (the concentration scales are different).

Simulated results of chloride concentration distribution in the WKP by the CFPv2 and UMT3D solute transport numerical model at the ends of the years 1967, 1986 and 2004 are discussed in this study. The concentration distributions at end of the model in 2018 as predictionare also included. At the end of 1967, the background chloride level of layer 1 in mainly part of the study site was about 2 mg/L and can be over 10 mg/L in northwest boundary of the study region. In layer 2, high concentration region at northwest and west boundary was much larger than layer 1, but chloride background level is lower than 1 mg/L in most of the region (Figure 6, left). Chloride concentration in the area close to SWF sprayfield significantly increased to higher than 5 mg/L in layer 1, which represented the effect of chloride leak from the SWF sprayfield facility and biosolids disposal site (Figure 6, right).

Figure 6.

Chloride concentration from the simulation by CFPv2 and UMT3D at the end of 1967 (left: layer 1; right: layer 2).

At the end of 1986, in addition to high chloride level inflow from the northwest boundary, the highest chloride concentration in layer 1 could be about 30 mg/L, had been observed at the SEF sprayfield (Figure 7, left). The chloride plumes with concentration higher 5mg/L can be found close to Lost Creek Sink, Wakulla Spring, Jump Creek Sink and Fisher Creek Sink. In layer 2, concentration near the west boundary could be as high as 15 mg/L, decreased to 6 mg/L near Wakulla Spring and lower than 3 mg/L at the east boundary (Figure 7, right). In the central of the study region, chloride plume area from the sprayfield also increased comparing with 1967 results.

Figure 7.

Chloride concentration from the simulation by CFPv2 and UMT3D at the end of 1986 (left: layer 1; right: layer 2).

Figure 8.

Chloride concentration from the simulation by CFPv2 and UMT3D at the end of 2004 (left: layer 1; right: layer 2).

At the end of 2004, simulated chloride level at the SEF sprayfield decreased to about 12 mg/L, which was much lower than the peak value (Figure 8, left). Concentration at the SWF sprayfield was about 5 mg/L due to the sprayfield shutdown. Chloride distribution in 2004 did not change much comparing with the 1986 distribution in the study region, indicated a relatively stable status. In layer 2, chloride plume at the SEF sprayfield was larger than that in 1986 and became an isolated high concentration zone (Figure 8, right). In other parts of the study region, chloride level also decreased from highest at the west boundary to lowest at east, which was similar to the distribution in 1986.

At the end of 2018 as a prediction of future, chloride levels in layer 1 at the SEF and SWF sprayfields keep decreasing since 2004 but are still about 10 mg/L and still higher than surrounding area (Figure 9, left). Chloride concentration at Spring Creek Springs is still the highest in the whole study region. In layer 2, chloride concentration in the north of the SWF sprayfield decreases below 5 mg/L because of the flushing effect (Figure 9, right). The west part of study region still has higher chloride concentration than the east part.

Figure 9.

Chloride concentration from the simulation by CFPv2 and UMT3D at the end of 2018 (left: layer 1; right: layer 2).

5. Discussion and conclusion

Groundwater flow in a well‐developed karst aquifer is mainly dominated within the fractures, conduits and caves formed and/or enlarged by carbonate dissolution. In general, most groundwater‐modeling methods, such as MODFLOW‐2000, apply Darcy equation to simulate groundwater flow, in which groundwater flow is laminar in the primary porosity (i.e., matrix porosity) of the aquifer. However, Darcian principle is not applicable in a well‐developed karst aquifer due to the dual permeability properties present in the aquifer. Therefore, advanced methodologies, such as the discrete‐continuum model, should be used to model groundwater flow in karst aquifers.

The physically based discrete‐continuum numerical models are widely applied to simulate groundwater flow and solute transport in a karst aquifer with conduit networks. The CFPM1 is able to accurately simulate groundwater flow in karst aquifers. The major hypothesis of this study was that CFPM1 would more accurately simulate flow in karst aquifers than MODFLOW‐2000 because CFPM1 accounts for turbulent flow, as well as pipe parameters such as pipe diameter, tortuosity, etc. In Gallegos et al. [22], a CFPM1 model of the Woodville Karst Plain was developed based on a previous MODFLOW‐2000 model by Davis et al. [35]. Both the CFPM1 and MODFLOW‐2000 models are created to run a transient simulation of an actual storm event. The simulated monitoring wells match well with the observations.

The coupled CFPv2 and UMT3D models are able to simulate solute transport in a karst aquifer. Chloride concentration plumes are originated from the sprayfields and extended southward to Wakulla Springs. Sprayfields in the south of Tallahassee and septic tanks in the rural area are the major source of chloride contamination in the WKP. From the monitoring wells observation, chloride concentration started to increased, especially after the SEF sprayfield went into service in 1980s and remained a high value in future prediction. Chloride contamination transported through the conduit network and formed extended plumes far away from the sprayfield. The transport method modified the flow and solute transport within conduits and considered the exchange between conduits and matrix, generated more accurate simulation results due to more precise description of conduits in the model.