**2. Study site and model setup**

nate aquifers supplyalmost 52%of allbedrockaquiferwithdrawalsonayearlybasis [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 ofthe retention andrelease effect between the conduits andmatrixdomains

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

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‐permea‐ bility 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],

[10, 11].

by the presence of a diffuse flow system.

116 Groundwater - Contaminant and Resource Management

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 MOD‐ FLOW‐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.

Numerical Simulation of Groundwater Flow and Solute Transport in a Karst Aquifer with Conduits http://dx.doi.org/10.5772/63766 119

conduits within the Woodville Karst Plain. The simulated heads, flow and transport within

In Gallegos et al. [22], the CFPM1 model of the Woodville Karst Plain is based on the MOD‐ FLOW‐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

the Woodville Karst Plain by Davis et al. [35] match with the observations very well.

118 Groundwater - Contaminant and Resource Management

numerical model: Wakulla Spring, Spring Creek Springs Group and St Marks Spring.

from onsite sewage disposal systems (e.g., septic tanks).

explored Wakulla caves system.

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

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

**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 con‐ ductivity 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.
