**3. Numerical simulation of groundwater flow in a karst aquifer**

#### **3.1. Governing equations of groundwater flow**

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

120 Groundwater - Contaminant and Resource Management

**Figure 3.** Pipe network in the CFPM1 subregional model of the Woodville Karst Plain. Modified from Gallegos et al.

[22].

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

$$\frac{\partial}{\partial \mathbf{x}} \left( K\_{\text{xx}} \frac{\partial h}{\partial \mathbf{x}} \right) + \frac{\partial}{\partial \mathbf{y}} \left( K\_{\text{yy}} \frac{\partial h}{\partial \mathbf{y}} \right) + \frac{\partial}{\partial \mathbf{z}} \left( K\_{\text{zz}} \frac{\partial h}{\partial \mathbf{z}} \right) \pm W = S\_s \left( \frac{\partial h}{\partial t} \right) \tag{1}$$

where *K*xx, Kyy and Kzz 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; *Ss* 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} \tag{2}
$$

where Δ*h* or *h <sup>L</sup>* is the head loss [L] measured along the pipe length Δ*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:

$$\mathcal{Q}\_{\alpha} = \alpha\_{j,\ i,k} \left( h\_n - h\_{j,i,k} \right) \tag{3}$$

where *Qex* is the volumetric flow exchange rate [L3 T-1]; *α <sup>j</sup>*, *<sup>i</sup>*, *<sup>k</sup>* is the pipe conductance at MODFLOW cell *j*, *i*, *k* [L2 T-1]; *h <sup>n</sup>* is the head [L] at pipe node *n* located at the center of the MODFLOW cell; and *h <sup>j</sup>*,*i*,*<sup>k</sup>* 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.

*xx yy zz s h hh h K K K WS x xy yz z t*

source term; *Ss* is the specific storage [L-1].

122 Groundwater - Contaminant and Resource Management

a conduit is conceptualized as a pipe,

gravitational acceleration constant [LT-2].

linear with head difference as follows:

MODFLOW cell *j*, *i*, *k* [L2

**3.2. Flow modeling application**

where *Qex* is the volumetric flow exchange rate [L3

where *K*xx, Kyy and Kzz 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/

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

> 2 *<sup>L</sup> l V hh f d g*

where Δ*h* or *h <sup>L</sup>* is the head loss [L] measured along the pipe length Δ*l* [L]; *f* is the friction factor [dimensionless]; *d* is the pipe diameter [L]; *V* is the mean velocity [LT-1]; *g* is the

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

> *Q hh ex j i k n j i k* = a

MODFLOW cell; and *h <sup>j</sup>*,*i*,*<sup>k</sup>* is the head [L] in the encompassing MODFLOW cell *j*, *i*, *k*.

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

2

<sup>D</sup> D= = (2)

, , ( , , ) (3)

T-1]; *h <sup>n</sup>* is the head [L] at pipe node *n* located at the center of the

T-1]; *α <sup>j</sup>*, *<sup>i</sup>*, *<sup>k</sup>* is the pipe conductance at

¶ ¶¶ ¶¶ ¶ ¶ æ ö æ ö æ ö æö ç ÷ + + ±= ç ÷ ç ÷ ç÷ ¶ ¶¶ ¶¶ ¶ ¶ è ø è ø è ø èø (1)

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.
