# 6.3 Distributed Parameter Models

Distributed parameter models may be deterministic or stochastic. For deterministic models, the output of the model is fully determined by the parameter values and the initial and boundary conditions. This is in contrast to stochastic models which have some type of inherent randomness (that is, the same input will lead to an ensemble of outputs because, for example, the input may specify statistics describing the frequency, size, and length of conduits, but not their exact size, length, or position, thus numerous realizations of the statistical distribution are generated to create individual deterministic simulations).

Various types of distributed parameter models have been applied to karst aquifers (Figure 69). The simplest approach is often referred to as the porous-equivalent media model approach (also called the single-continuum porous-equivalent [SCPE], heterogeneous-continuum, distributed-parameter, or smeared-conduit model). The SCPE assumes that flow in the aquifer system can be simulated with the potential flow equation of fluid mechanics which is the way typical porous aquifers are represented.

Figure 69  Approaches to karst modeling with distributed parameter models (numerical models) ranging from porous-media equivalents to fracture and conduit networks that are simulated with flow equations for fractures or conduits. Stochastic methods are commonly used for defining the nature and location of fractures and conduits. Modified from Teutsch and Sauter (1998).

Regardless of the distributed parameter modeling approach, field data are collected to: map the aquifers and confining units; determine the location and quantity of water entering and exiting the system; and, provide water-level and flow observations for calibration of the model. The data are used to define the hydrogeologic framework and develop a conceptual flow model for the application and finally for calibration of the site model. Calibration is the process of adjusting the value of model parameters until the values simulated by the model match the measured field observations. It is best if a basic model is developed early in the investigation because the process of building a model can reveal the type and location of data that should be collected during the rest of the investigation in order to improve the model. A common mistake is beginning the modeling task after the field work has been completed, at which time it is too late to use the model to inform the investigation process.

Most distributed parameter models have limitations in their use and application resulting from: simplifications of the system, inadequate calibration data, or poorly understood system geometry and boundary conditions. It is beyond the scope of this book to discuss the documentation of models such that their limitations are fully presented. However, good resources on this topic can be found in Anderson and others (2015) as well as in Reilly and Harbaugh (2004).

Distributed parameter models are frequently applied to karst aquifer systems (Kuniansky and Holligan, 1994; Teutsch and Sauter, 1998; Kuniansky et al., 2001; Kuniansky and Ardis, 2004; Scanlon et al., 2003). These types of models are process oriented, and the process can be represented by the Navier-Stokes equations that describe unsteady fluid flow in three dimensions. The Navier-Stokes equations are coupled sets of partial differential equations that describe how pressure, temperature, density and viscosity of a moving fluid are related. The equations were derived based on the underlying principles of conservation of mass, energy, and momentum and are the basis for almost all fluid mechanics problems (Daily and Harleman, 1966).

Unfortunately, all the forces described by these partial differential equations cannot be solved readily on the computers that are available at the time this book is being written due to speed, memory and storage limitations. Consequently, forces and processes that have negligible impact on the flow problem can be eliminated from the full set of Navier-Stokes equations such that a reduced equation (or set of equations) that includes the important physical process(es) is (are) applied. In groundwater problems, flow is typically dominated by potential energy as reflected by the hydraulic head gradient within an aquifer. The energy that would go into eddies or vortices within the water (kinetic energy losses from turbulence or non-laminar flow) is negligible for most groundwater flow problems as are differences in temperature, density, and viscosity for shallow groundwater systems. Thus, the potential-flow equation, used in groundwater model codes such as the basic version of the United States Geological Survey software called MODFLOW (McDonald and Harbaugh, 1988), assumes that flow is laminar and that temperature, density, and viscosity are constant over the model domain. Provided that the problem is reasonably well represented by these assumptions, the computational burden and time required to run the simulation can be greatly reduced with little loss of accuracy in the solution.

The potential-flow equation has also been used in aerodynamics for flow of air over a wing; in thermodynamics for heat flow; in hydrodynamics for flow around objects in a stream; for full pipe flow; and for flow of electrical current. Distributed parameter models require dividing (discretizing) the aquifer system into representative three-dimensional volumes, such as finite-difference cells or finite elements. Depending on the scale of the study, representation of the karst aquifer may also be greatly simplified due to limitations on the total number of cells or elements that can be solved in a reasonable time on a computer. Thus, the greater the area and thickness of the simulated karst system, the larger the number of cells or elements; which may be less representative of fine-scale features and require the use of composite hydraulic properties (Sepúlveda and Kuniansky, 2009). The physical aquifer properties (parameters), such as storage and hydraulic conductivity are distributed spatially and may have a different value in each model cell or element. Thus, the general name, distributed parameter model, is widely applied to these mathematical modeling approaches.

Other distributed parameter approaches shown in Figure 69 are the dual-continuum porous-equivalent model (DCPE), the hybrid model (HM), the discrete single fracture or conduit set model (DSFS), and the discrete multiple fracture or conduit set model (DMFS). The DCPE model links two flow regimes (two SCPEs) at each cell or element with a head-dependent water exchange term between sub-cells or sub-elements. In a DCPE the flow within each cell represents a physical volume of the aquifer with two, sub-cells exchanging flow with one another. One sub-cell represents the conduits and the other sub-cell represents flow in the rock matrix. The hybrid model (HM) couples a three-dimensional SCPE model with a discrete one-dimensional flowing fracture or conduit network (Figure 69, HM). The fracture-set or conduit-network modeling approaches shown in Figure 69 may involve the stochastic generation of fracture or conduit networks that are simulated with equations for conduit or fracture flow. The discrete, single-fracture or conduit network set (DSFS) and discrete multiple-fracture or multiple-conduit network set (DMFS) models are rarely applied to karst aquifers. The advantages and disadvantages of each modeling approach are described in the following sections.

Most distributed parameter models of karst systems have a large number of cells or elements used to subdivide the aquifers and confining units into discrete cell volumes (discretization). For illustrative purposes, Figure 71 shows the nomenclature for the division (discretization) of a five-layer aquifer system into cell volumes used by the block-centered finite-difference groundwater flow modeling code MODFLOW. Figure 71 is a simplified model diagram where all layers are constant thickness, and the number of rows and columns is small. Most site models have hundreds of rows and columns and the row and column dimensions are not constant.

Figure 70  Distributed parameter models divide (discretize) the aquifer system into representative three-dimensional volumes and define properties that represent the equivalent bulk properties of the volume of rock that the block represents. This diagram shows the convention for numbering block-centered, finite-difference cells for MODFLOW. In the illustrated case, each vertical layer has a different lithology as indicated by the shading pattern. Modified from McDonald and Harbaugh (1988).

# Single-Continuum Porous-Equivalent Models

An SCPE approach using the potential-flow equation (which assumes laminar flow at constant temperature, density, and viscosity) is the simplest distributed parameter model to apply. Figure 72 displays a simple example of applying a SCPE to karst using MODFLOW. The model has two layers. The top layer is predominantly unconsolidated sand, and the bottom layer is a confined karst aquifer with mapped conduits that connect to an artesian spring that discharges in a reach of a perennial stream in the top layer. It is simplified by using high hydraulic conductivity cells to represent conduits.

The SCPE approach has been applied for regional, or sub-regional, flow in some karst aquifers because the scale of investigation is much greater than the scale of the heterogeneities in the rocks. In this case, heterogeneities refer to dissolution features or flow units with different water-transmitting properties. The SCPE approach has been successfully applied for water-resources investigations where the models are calibrated to steady-state average conditions (generally representing one season- or annual-average recharge and pumpage conditions during the time period when the potentiometric mapping was conducted). The SCPE approach has also been applied for transient conditions using multiple stress periods that may represent monthly, seasonal, or annual flow. The models simulate head (water level) and provide a regional- or sub-regional-scale water budgets for water planning purposes. In general, the SCPE approach can simulate transient average monthly or annual spring flow, but generally cannot reproduce detailed storm-event hydrographs as well as other model types because non-laminar flow occurs during storms (Hill et al., 2010; Kuniansky et al., 2011; Gallegos et al., 2013; Saller et al., 2013; Kuniansky, 2014, 2016).

Simulations of advective transport using a single-continuum model are infrequently performed and have varying degrees of success. Advective transport modeling within karst aquifers has been able to match estimated geochemical age and travel times derived from tracer-tests for the Floridan Aquifer Systems in the USA by using effective porosity of less than 5 percent rather than the total porosity which ranges from 15 to 40 percent. These models have been documented by Knochenmus and Robinson (1996) Kuniansky and others (2001), Merritt (2004), Renken and others (2005), as well as Davis and others (2010). In these studies, groundwater flow volume (Darcy flux) was mainly governed by the model cell hydraulic conductivity and model layer thickness. Effective porosity is used for calculation of a pore velocity that matches model simulated values to observed water age or time of travel determined by using field tracer tests. In some cases where conduit locations were known, the finite-difference cells or finite elements with conduits were assigned much greater hydraulic conductivity values than surrounding cells and were successful in reproducing transient spring discharge (annual, seasonal, and/or monthly averages) and matching tracer-test-derived times of travel (Davis et al., 2010). A similar attempt was not as successful within the Edwards Aquifer system in Texas, USA (Lindgren et al., 2004, 2009), but was successful in the Fort Payne aquifer in Tennessee, USA, which is a less permeable karst aquifer than the Edwards Aquifer system or the Florida Aquifer System (Haugh, 2006).

Figure 71  Example of a single-continuum porous equivalent model. This is a two-layer model with karst beneath a uniform, 100 m thick sand. All cells are 100 by 100 m on each side. Limestones are represented in layer 2 with cell hydraulic conductivity ranging from 100 to 10,000 m/d. Layer 2 is 200 m thick. Mapped conduits are located in cells with K of 1,000 and 10,000 m/d. The mapped conduit cells, shown in red and yellow, are mostly in layer 2, with one cell in layer 1 where the major spring outlet occurs and connects to the cell below in layer 2. A perennial stream flows across layer 1 with a stage of 95 m at row 1 column 30 and 94 m at row 30 column 1. The river cells are indicated with a light blue line within each cell in layer one and the river stage decreases linearly downstream. A constant rate of recharge is applied to the top of the model.

# Dual-Continuum Porous-Equivalent Models

The DCPE models link two, single-continuum, porous-equivalent models (SCPE) via a head-dependent flux term between the linked cells of the SCPEs. Thus, flow in each aquifer cell volume is represented by two, sub-cell properties of each SCPE within the dual porosity karst system, one represents flow in conduits and the other flow in the rock matrix. DCPE models were designed to help subdivide the flow within a block of aquifer between the conduits and the rock matrix by assigning different hydraulic properties to each, such that the faster flow in the extremely transmissive conduits and the slower flow in the rock matrix could result in better matches to observed spring flows.

In simple terms, a head-dependent flux means that flow between the two, sub-cells in the dual-continuum cell is controlled by the head difference between the two continuums and the hydraulic conductivity between them (water flows from the continuum with high head to the continuum with low head). One SCPE has large hydraulic conductivity and small storage, representing conduits, and the other SCPE has smaller hydraulic conductivity and larger storage, representing the rock matrix. The groundwater flow equations are solved iteratively for head at all SCPE sub-cells of the DCPE model for each step that the model takes through time until the solution converges. Most DCPE software codes look similar to a block-centered, finite-difference, numerical model for SCPE (Figures 70 and 71), but with two arrays of rock properties in the layer(s) that are defined as having dual porosity.

DCPE models have been applied successfully when the geometry of conduits is not known (Teutsch and Sauter, 1991; Sauter, 1993; Teutsch, 1993; Teutsch and Sauter, 1998; Painter et al., 2007). One advantage of using the DCPE approach is that the detailed geometry of the conduits is not required (Teutsch, 1993; Sauter, 1993; Lang, 1995). The main advantages of using a DCPE model are its capacity to simulate rapid variations in discharge and head following recharge events and its mathematical representation of both the conduit and rock matrix contribution to total flow through time. The required data for the model is relatively modest, and the effort required to construct the model is manageable for most projects. The DCPE model, however, generally does not have the ability to simulate transport processes on a small scale (Mohrlok, 1996). Application of this modeling approach is not common.

# Hybrid Models

A hybrid model (HM) is the coupling of an SCPE model with a discrete one-dimensional conduit network model (Kiraly, 1998; Teutsch and Sauter, 1998). The HM approach allows the integration of detailed information about conduits in areas where the geometry of major conduits is known, thus providing a more representative model of the physical system. The HM links a three-dimensional SCPE model for groundwater flow within the rock matrix with a one-dimensional conduit network model that can have laminar or turbulent flow. The conduit network geometry is specified by defining locations within finite-difference cells where conduits are connected to other conduits (conduit network nodes). The SCPE model exchanges water between the conduit network and the porous matrix in cells containing conduit nodes by using head-dependent flux terms and iterates calculations between the two models until both converge to a solution for head, in a manner similar to the DCPE model.

The 2008 version of the MODFLOW-2005 conduit flow process (MODFLOW-CFP) (Harbaugh, 2005; Shoemaker et al., 2008) is based on the hybrid modeling approach developed in the Carbonate Aquifer Void Evolution (CAVE) software code (Clemens et al., 1996; Bauer et al., 2000; Bauer, 2002; Birk, 2002; Bauer et al., 2003; Liedl et al., 2003). The main advantage of this modeling approach is that it can simulate high transport velocities, under laminar and non-laminar flow, as observed in karst aquifer system conduits; while accounting for the presence of a lower hydraulic conductivity rock matrix. Most of the water storage occurs within the karst aquifer rock matrix. This approach has limited application for most studies due to the lack of geometric data for the entire conduit system (conduit location, diameter, roughness, and tortuosity) given that often only the largest conduits are mapped. The first release of the CFP has limitations in that the conduit network is intended for conduit-full flow where flow is controlled by the pressure gradient (potential flow but including turbulence) and does not account for free-surface flow in partially full conduits. For free-surface, open-channel flow, the free water surface is in equilibrium with atmospheric pressure and the forces driving flow are dominated by gravity. In this case, the gravity terms cannot be omitted from the Navier-Stokes equations because flow is predominantly controlled by open-channel slope, roughness, and channel shape. In short, the 2008 version of CFP should not be used for systems with large networks of partially full conduits. MODFLOW-CFP and research variants of the code have been applied to the Florida Aquifer System in several locations as discussed in the example for this section with respect to an application in the Woodville Karst Plain, Florida, USA.

The MODFLOW-CFP simulation code includes three modes. Mode 1 is the hybrid model described above. Mode 2 includes the capability to insert a high-conductivity flow layer that can switch between laminar and turbulent flow. Mode 3 allows for coupling a discrete one-dimensional conduit network in a high-conductivity flow layer that can switch between laminar and turbulent flow. Figure 72 shows how one-dimensional pipe conduit networks can be configured within the MODFLOW cells. A maximum of six pipes can be connected at a single node. Even if the mapped conduit is not at the centroid of a cell volume, the head-dependent flux term between the MODFLOW cell and the pipes is calculated between the pipe node and centroid of the cell. If the modeler wants to better locate the connection, then a finer grid needs to be defined. Additionally, the pipe leakage conductance is calculated using half of the surface area of the pipe between two nodes. Shoemaker and others (2008) provide more detailed discussion of MODFLOW-CFP.

Figure 72  Three examples of how the one-dimensional conduit network could link to adjacent finite difference cells in MODFLOW-CFP. A maximum of six pipes can connect at a single node and flow exchange between the pipe and the model occurs at nodes within the specified finite-difference cell. Modified from Shoemaker and others (2008).

Models originally intended for surface-water systems have been used to represent karst aquifers. The Storm Water Management Model (SWMM) (Metcalf and Eddy Incorporated, 1971) was designed for simulation of sewer systems and was applied by Peterson and Wicks (2006) to simulate karst conduits, but no interchange with the rock matrix was simulated. The MODBRANCH code (Swain and Wexler, 1996) was modified for simulation of a karst system by Zhang and Lerner (2000). MODBRANCH was originally developed by coupling one-dimensional, free-surface, open-channel flow with the top layer of a three-dimensional groundwater flow model for simulation of interaction between groundwater and surface-water. Reimann and others (2011b) developed the ModBraC version of a hybrid model that couples the SCPE model with a conduit network model capable of simulating storage of the conduit system, conduit-full flow, and open-channel flow in the conduits. Accurate representation of conduit storage is essential in order to simulate the lag time between discharge changes and indicators of transport (for example, temperature, conductivity) after recharge events. Grubbs and Crandall (2007) applied MODBRANCH for simulation of groundwater and surface-water interactions in a karst aquifer.

# Discrete Single and Multiple Fracture or Conduit Network Models

The discrete fracture approach using a single fracture set (DSFS), or multiple fracture sets (DMFS) has been proposed for problems where transient solute-transport responses are desired for an aquifer system dominated by fractures or conduits (Adams and Parkin, 2002). Knowledge of the fracture network geom­etry is required for application of these models. Generally, data from mapping of the fracture network(s) are limited, thus stochastic methods are typically used to generate the single- or multiple-fracture network for individual, deterministic, numerical simulations of flow. The disadvantages of discrete multiple fracture networks are the requirement of detailed knowledge of a fracture network at multiple scales and the application of computationally intensive codes with long computer simulation time (days) and large memory requirements (Lang, 1995). Much effort is in the collection of data on fractures for stochastic generation of multiple fracture networks from the statistics based on fracture data from borehole imagery and processing of the images to estimate aperture size and orientation. Often fractures are planar features within low-permeability hard rock and borehole image processing is one of the best tools to delineate their orientation and gather the statistics for stochastic model simulations. Many different algorithms are used to generate the fractures and convert these to bulk cell-block properties. Lei and others (2017) provide an excellent overview of discrete fracture network modeling. Implementation of fluid-flow simulation is usually similar to that of SCPE or DCPE with rock properties assigned to cells based on the fracture network modeling (Figure 73); however other methods are applied in some cases. Application of this approach to field-scale problems is not common, so no examples are available. The petroleum industry has applied DSFS/DMFS, but model run times have not been published (McClure and Horne, 2013).

Figure 73  Examples of fracture network modeling. a) A network is stochastically generated for each unit, then a grid is overlain on that layer and bulk properties are calculated for each finite difference cell for dual or single continuum flow simulation. b) Fractures are mapped (those oriented vertically are in red and those oriented horizontally in green) and the cells of the finite-difference grid are assigned properties based on the mapped fractures that intersect most of the finite-difference grid. Modified from Lei and others (2017).

# Example Comparisons of Single-Continuum and Hybrid Models

Extensive mapping of the sinks and submerged conduit system of the Wakulla Springs-Leon Sinks, Florida, USA, has been accomplished over the past 20 years by cave divers of the Global Underwater Explorers as part of the Woodville Karst Plain Project (WKPP) as described by Kernagis and others (2008). As a result of their pioneering work, this karst drainage system is believed to be well characterized. The Woodville Karst Plain is 450 square miles (~1200 km2) and contains the Wakulla Springs-Leon Sinks system. Wakulla Springs is one of the largest first-magnitude springs in Florida, with average flows of approximately 400 ft3/s (cubic feet per second; ~11 m3/s). The spring occurs at a large vent opening of about 25 m by 15 m which forms a pool that is approximately 90 m in diameter and 55 m deep at the headwaters of the Wakulla River at Edward Ball Wakulla Springs State Park, Florida, USA (Florida Department of Environmental Protection, 2016). A satellite image of Wakulla Springs is shown in Figure 1c. This area has been investigated for decades and several SCPE models have been developed and calibrated.

The Woodville Karst Plain Project (WKPP) maps were used by Davis and others (2010) to update an older SCPE model that had been developed and calibrated for a previous study (Davis and Katz, 2007). A hybrid model was developed and calibrated based on the model developed by Davis and others (2010) and published by Gallegos and others (2013). Ultimately, three modeling approaches were developed and their results were compared to examine the applicability of the approaches for this area.

Background on the Three Modeling Approaches

Davis and others (2010) used a SCPE approach to simulate transient conditions in the area shown by the sub-regional model boundary in Figure 74. The simulation begins on January 1, 1966, when spray-field operations (irrigation using treated wastewater as part of municipal wastewater treatment processing) began, and ends in 2018, when operating system upgrades to the wastewater system were planned to take effect. Changing hydrologic stresses (for example, pumping and recharge) on the groundwater system were primarily represented in the model as average-annual conditions throughout the 1966–2018 model simulation period. For a 2.2-year period during 2006 and 2007 when tracer tests were conducted by Hazlett-Kincaid, Incorporated, hydrologic conditions were represented by average conditions over 10-day periods. Data describing the hydrologic stresses are listed in Davis and others (2010). Observations of water levels and spring discharge were available for November 1991 and May to early June 2006. These data were used to calibrate the model.

Figure 74  The extent of the sub-regional model of the Woodville Karst Plain, Florida, USA with colors indicating the hydraulic conductivity in model layer 2. Large hydraulic conductivities correspond to submerged conduits. Modified from Davis and others (2010).

Two additional models were developed such that three approaches were used to represent the system. One, herein called Approach 2, was an SCPE model as used by Davis and others (2010), but with both laminar and turbulent flow in the SCPE for layer 2. The other, herein called Approach 3, was a hybrid model (HM) consisting of a single continuum as used in Approach 1 and 2 coupled to a one-dimensional, conduit-flow network capable of simulating laminar and non-laminar flow in the conduit network.

The three modeling approaches are summarized here to help the reader keep their differences in mind.

• Approach 1—SCPE model with only laminar flow using a previously calibrated model of Davis and others (2010). This model has high hydraulic conductivity (K) cells in mapped and inferred conduits.
• Approach 2—SCPE model as in Approach 1, but with both laminar and turbulent flow in the SCPE for layer 2. This is an application of MODFLOW-CFP using mode 2 (Shoemaker et al., 2008; Kuniansky et al., 2008; Davis et al., 2010; Kuniansky et al., 2011; Reimann et al., 2011a, b; Kuniansky, 2014, 2016). Groundwater flow can be laminar or non-laminar within a SCPE model layer for the MODFLOW-CFP mode-2 approach. Water temperature, average pore diameter, and critical Reynolds numbers for layer 2 are parameters that were not needed for Approach 1, but were required for Approach 2.
• Approach 3—HM model consisting of a single continuum coupled to a one-dimensional, conduit-flow network capable of simulating laminar and non-laminar flow in the conduit network. This is an application of MODFLOW-CFP using mode 1 (Shoemaker et al., 2008; Gallegos et al., 2013; Kuniansky, 2014, 2016). This model was developed from the Davis and others (2010) model with the hydraulic conductivity in their high-K cells decreased to the lower background-K values of surrounding cells and an interconnected pipe network was used to represent conduits.

More details about the hydrogeology of the area as well as the original model grid, boundary conditions, and calibration data, which were used as the basis for calibrating all three approaches, are discussed in Davis (1996), Davis and Katz (2007), Davis and others (2010), as well as Gallegos and others (2013). The simulation input and output files discussed herein are available as a data release (Kuniansky, 2016b).

The models for all three approaches have the same two-layer discretization (Figure 75). The top layer is the confined upper 60 m of the Upper Floridan aquifer, and the second layer is the confined lower part of the Upper Floridan aquifer, which is more than 300 m thick in the west, thins to 60 m in the east, and contains the submerged mapped conduits (Davis et al., 2010).

Figure 75  Generalized geologic cross section and model layers for the Woodville Karst Plain models, Florida, USA. From Davis and others (2010).

The calibrated K values used by Davis and others (2010) are shown in Figure 74, with hydraulic conductivity values greater than 10,000 ft/day (feet per day; > 3000 m/d) assigned to cells that represent mapped and hypothesized conduits. The same values were used for SCPE Approach 2. The only additional parameter values required for Approach 2 were an average pore diameter of 0.75 ft (0.23 m), and lower and upper critical Reynolds numbers of 11 and 150, respectively. After Davis and others completed their project, head and spring flow data were collected for a 52-day storm event in 2008, so the pore diameter and critical Reynolds numbers were adjusted in a transient calibration of the SCPE Approach 2 model using average-daily hydrologic conditions. For model Approach 2, the average pore diameter and the upper and lower critical Reynolds number values were adjusted to improve the match of simulated and observed flows at Wakulla Springs for the daily transient simulation of the 52-day storm. The MODFLOW-CFP used for Approach 2 allows non-laminar flow, but is simplified by requiring only temperature, average pore diameter, and upper and lower Reynolds numbers for the entire layer so adjustments to improve the match at one spring impacts the simulated flow at other springs as well, thus the calibration strove to only match discharge at Wakulla Springs. Temperature was held constant at the average groundwater temperature (20.5 oC). In this original version of CFP mode 2, the effect of non-laminar flow requires specification of parameters that represent large-pore porous media, not conduits. Reimann and others (2011a) developed a modification that allows for specification of Reynolds numbers more typical of large conduits. When the Reimann and others (2011a) form of CFP mode 2 was applied to this system, it reproduced similar spring flow as the calibrated model of Approach 2 by using higher values (an average pore diameter of 1 foot (~0.3 m), and lower and upper critical Reynolds numbers of 3,000 and 10,000, respectively).

Model Approach 3 is the HM with the one-dimensional, conduit-flow network representing the conduits. It was developed by replacing the higher K values of the SCPE model (those greater than 3,000 m/d) with the smaller K values of the surrounding cells and adding the mapped conduit network to those cells as one-dimensional circular pipe flow. Unlike the Approach 2 model, the Approach 3 model had to be calibrated because the initial parameter values used to describe the conduit properties did not produce an acceptable match to the calibration data. During the calibration process, the background K’s were not changed from those used in Approach 1 and 2, only the values of the pipe parameters were adjusted to obtain an acceptable fit to the calibration data.

The hybrid model is large with over 1000 interconnected pipes and nodes. The time required to run the HM transient model is ~100 times longer than the SCPE models. This makes calibration extremely difficult because the model needs to be run every time a parameter value is adjusted, and values need to be adjusted many times to obtain an acceptable fit to the field observations. Consequently, this model was not calibrated to the transient data used by Davis and others (2010). It was only calibrated to steady-state scenarios that approximated average conditions for two periods when both spring discharge and potentiometric data were available including:

• November 1991
• May-June 2006

Further calibration of pipe properties was accomplished by matching the times of travel along conduits from dye tracing tests conducted during 2006 using the steady-state data for May-June 2006 from the following sink-to-spring dye-trace times. Davis and others (2010) and Gallegos and others (2013) provide more details as both models were calibrated to match these dye-trace travel times:

• Fisher Sink to Wakulla Springs took 10 days;
• Ames Sink to Wakulla Springs took 20 days; and,
• Turf Sink to Wakulla Springs took 40 days.

Even steady-state HM simulation requires many iterations before the SCPE and the one-dimensional pipe network reach what is called convergence where the simulated heads and flows are not changing with additional iterations. Gallegos and others (2013) provide details about the pipe-network parameters and the calibration. Only the parameter values describing the pipe network were adjusted when calibrating the HM model.

The modelers defined criteria for the calibration, requiring that parameter values be adjusted until:

• simulated spring flows were within 10 percent of measured flows; and,
• simulated heads were within plus or minus 2 m of measured heads.

The SPCE models (Approaches 1 and 2) and the HM model (Approach 3) all simulated observed average spring discharge at Wakulla and Spring Creek Springs to be within the calibration criterion of 10 percent. For all three approaches, the difference between simulated heads and heads measured in the field were within +/-2 m (Davis et al., 2010; Gallegos, 2011; Gallegos et al., 2013). For average spring discharge, the three modeling approaches produce similar results, and all are considered acceptable representations of the system (Davis et al., 2010; Kuniansky et al., 2011; Gallegos et al., 2013; Kuniansky, 2014, 2016).

While all three modeling approaches met the water-level and flow calibration criteria, the simulated potentiometric maps for model layer 2 (the layer with submerged conduits) differed slightly between the HM (Approach 3) and the two SCPE models (Approach 1 and 2) in that, although the simulated heads were within +/-2 m of the measured heads, the simulated heads in layer 2 of the HM were nearly all lower than the measured heads. That is, the head differences were biased (Gallegos et al., 2013) thus, it was a less-than-acceptable representation of the system. Some decrease of all the hydraulic conductivities or perhaps only those of the matrix or the conduits is likely to resolve the bias, but that was not pursued in the study. The calibrated parameters would be different if the HM was rigorously calibrated to match the larger transient data sets rather than the steady-state average conditions, but transient calibration was not possible with the available computing power.

Scenarios Simulated with the Calibrated Models

Three scenarios were simulated and compared using the various model approaches to examine applicability of the approaches for representing the Floridan Aquifer System in the vicinity of Wakulla Springs. The three scenarios represent the following conditions.

• Scenario 1 – transient representation of spray-field operations as presented by Davis and others (2010) for the 2.7-year period with 10-day averaging of hydrologic conditions (81 10-day average transient conditions beginning January 1, 1966 and extending 2.7 years).
• Scenario 2 – steady-state representation of average conditions for the entire period from 1966 to 2018.
• Scenario 3 – transient representation of a 52-day storm that occurred in 2008 using average-daily hydrologic conditions. This was a storm event that started on August 13, 2008 for which daily discharge was gaged at Wakulla Spring and daily recharge and groundwater pumpage was estimated for the simulations.

Comparison of transient representation of spray-field operations

Scenario 1 used the calibrated model for Approach 1 (SCPE with high K cells in layer 2 where conduits are mapped or inferred) fully documented in Davis and others (2010). Scenario 1 was conducted with Approach 2 only (SCPE with turbulence allowed in layer 2 and same hydraulic conductivity and storage properties as Approach 1). Small differences occurred between heads and flows when turbulence was turned on with the average conditions used by Davis and others (2010). With annual and 10-day averaged recharge and pumpage, turbulence occurred in the high-K model cells (cells with conduits) but did not change model head results enough to be noticed on potentiometric maps even with 5 ft (1.5 m) contour maps. For almost all the 10-day periods turbulence did not extend beyond the cells with conduits. For two of the 10-day periods with larger recharge (periods 44 and 66), turbulence occurred in the more cells than any other period and in a few adjacent to the cells with conduits. While it would be difficult to notice any difference in simulated water levels, there was a minor difference in the cumulative water budget from the comparison of Approach 1 and 2. The small amount of turbulence resulted in slightly less flow through the system for the 2.7-year period. The reduction was 8 percent, which is negligible given that the estimates of recharge and pumpage have more uncertainty than 8 percent.

Comparison of Simulated Potentiometric Maps for Steady-State Average Conditions

For illustration purposes, three simulated potentiometric maps are shown for model Layer 2, Scenario 2, which is a steady-state condition created by averaging the hydrologic stresses of recharge and groundwater withdrawals from 1966 to 2018 (Figure 76). The potentiometric map in Figure 76a shows a poor representation of the system with homogeneous K in layer 2 and no high hydraulic conductivity cells where conduits are mapped in an SCPE model. Figure 76b shows the essentially identical head distributions obtained from Approaches 1 and 2 which are SCPE models with high hydraulic conductivity cells where conduits are mapped (Figure 74 shows the high K cells with conduits). Figure 76c shows the simulated potentiometric map from the HM (Approach 3). The simulated potentiometric surfaces are similar when the conduits are represented in either an SCPE model or a hybrid model (Figure 76b and c). Thus, incorporation of conduits with very large hydraulic conductivity cells in an SCPE approach can mimic the HM approach. If high K cells are not used in areas of conduits, then the SCPE model will not approximate the water levels of an HM model. Turbulent conditions were not a factor for the average steady-state conditions. In general, scuba divers explore and map the submerged conduit system during low or base flows. Even though the low flow rates at Wakulla Spring are large, often over 100 ft3/s (3 m3/s), scuba divers do not enter these springs during storms and they have not observed turbulent conditions at these low or average flows. This can be observed in the video link associated with Figure 39 that shows scuba divers in the Weeki Wachee Spring conduit system.

Figure 76  Steady-state potentiometric maps of model layer 2 for simulation of average conditions from 1966-2018 in the Woodville Karst Plain, Florida, USA. a)  Simulated potentiometric surface for an uncalibrated single-continuum porous-equivalent media model with no high hydraulic conductivity cells to represent conduits. b) Simulated potentiometric surface for the calibrated single-continuum porous-equivalent media model of Approach 1 with higher hydraulic conductivity in cells thought to contain conduits as indicated by small blue rectangles. There was no visible difference in heads between model Approaches 1 and 2 because, although turbulence could occur when using Approach 2, it did not occur for the simulated conditions. c) Simulated potentiometric surface for the calibrated hybrid model (model Approach 3) with low hydraulic conductivity values representing the rock matrix and a one-dimensional conduit network in cells with conduits. The equivalent porous media approaches (SCPE) produced a head distribution that was nearly identical to the distribution produced by the hybrid model approach (HM) as indicated by the similarity of contours in b) and c). Modified from Kuniansky (2016 a,b).

Comparison of Transient Simulations of a Storm Event using AverageDaily Conditions

In order to observe larger differences between the SCPE and hybrid model approaches, a transient period with daily spring flow observations and a storm event (rising and falling of spring flow) was simulated without recalibrating the models to the storm event data (Kuniansky et al., 2011; Gallegos et al., 2013; Kuniansky, 2014). The 52-day period began August 13, 2008 and was simulated with daily stress conditions. Spring discharge rose rapidly on August 22 (10th day) as shown in Figure 77 and peaked on August 26 (14th day) at 2,363 ft3/s (cubic feet per second; ~67 m3/s) as indicated on Table 8. Only the Approach-2 model simulated the observed peak discharge within 10 percent of the measured discharge (Figure 77). The shape of the storm hydrograph for Wakulla Springs is best matched by model Approach 2, the SCPE model with turbulence (Figure 77).

Figure 77  Simulated and observed spring flow at Wakulla Springs, Florida, USA. Unit conversion: 1000 ft3/s (cubic feet per second; ~28.3 m3/s). Modified from Kuniansky (2014).

Table 8  Peak discharge at Wakulla Springs, August 26, 2008 during a storm event

 cubic feet per second cubic meters per second percent difference from observed comment Observed spring discharge 2,363 67 Approach 1, SCPE no turbulence 3,085 87 -31 overshoots peak Approach 2, SCPE allowing turbulence 2,149 61 9 reasonable match Approach 3, Hybrid Model 1,653 47 30 undershoots peak

Model Approach 1 overestimated and model Approach 3 underestimated the peak daily discharge at Wakulla Spring and both were within 30 percent of the observed peak value. Model Approach 2 underestimated the peak but was within 9 percent of the observed peak daily discharge. When no turbulence is considered in the modeling approach then discharge remains a linear function of head gradient and thus Approach 1 overshoots the peak discharge. When turbulence is incorporated into the models, then discharge does not increase linearly with the head gradient and so the simulated peak discharge is not as high. In both Approach 2 and 3 the peak discharge is reduced by the turbulence. With a storm event and daily hydrologic stresses applied to the simulations, turbulence is a factor in both Approach 2 and 3 and reduces the peak discharge.

When the total volumes of measured and simulated discharges for the 52-day period at Wakulla Springs were compared, the following was found:

• Approach 1 (SCPE laminar) matched observed volume within 23 percent;
• Approach 2 (SCPE laminar and non-laminar) matched within 17 percent; and,
• Approach 3 (HM) matched within 0.01 percent.

Simulated fits to observed flow at Spring Creek Spring and Saint Marks Spring were worse than those at Wakulla Springs for all three modeling approaches and those results are not shown herein. Given that none of the approaches captured both the peak and the total volume, none of these models are satisfactory for simulation of this storm event. All require more calibration. When considering this, it is useful to remember that, with the exception of the average pore diameter and critical Reynolds numbers of Approach 2, the models were not calibrated using average-daily stresses to match the daily discharge data. Models of Approach 1 and 2 were calibrated to transient hydrologic conditions, but these were mostly annual-average conditions with only about 2 years of average 10-day conditions. The model of Approach 3 was only calibrated to steady state conditions because of computational limitations. Calibration to daily discharge data using average-daily stresses would improve predictive capability of the models.

The computation time required by the HM Approach 3 to simulate the 52-day transient period (about 12 hours in 2016) was more than 100 times longer than the time required for the SCPE model Approaches 1 and 2. Additionally, the run times for the hybrid model were long and thus it was not calibrated for the full transient conditions from 1966 to 2018 as used for model Approaches 1 and 2. This is a large hybrid model, with more than 1,000 simulated conduits and nodes. The hybrid approach requires that two models are solved iteratively until both converge for each 1-day step through time, requiring a lot of computation time. Practical model applications need to balance the rigor of the solution and the accuracy of the predictions.

Other studies have shown more success with HM. For example, Hill and others (2010) compared a hybrid model approach with a SCPE model approach for two other spring systems in Florida (Weeki Wachee and Twin Dees) and concluded that the hybrid model more closely simulated observed transient spring discharge. Another study (Saller et al., 2013) converted a SCPE model to a hybrid model for a watershed in the Madison aquifer in South Dakota and found that the hybrid model better matched head observations at monitor wells, and both models simulated spring discharge within the calibration criteria.

Although other studies have found HM to provide better results, the Wakulla Springs example indicates that, for simulation of average conditions, none of the three approaches is distinctly better. However, the HM was manually calibrated with steady-state periods rather than the full transient simulation from 1966 to 2018 that was used for the SCPE models owing to long run times, which put the approach at a disadvantage (Gallegos et al., 2013). The HM is likely to have provided better results if it had been calibrated to the same transient conditions.

The SCPE model with turbulence (Approach 2) best matched peak spring daily discharge at Wakulla Springs for the 52-day storm, and the HM (Approach 3) best matched the total volume of spring flow. None of the model approaches had acceptable matches to the observed daily storm hydrographs at the springs in the model domain. The model of Davis and others (2010) was not intended for matching storm events. This illustrates the important concept that, if a project requires accurate simulation of daily storm events, then the model must be calibrated using storm event data.

It is unlikely that the extra effort required to use a hybrid model, both in data preparation and computation time, is justified for investigations not requiring site-scale transport or prediction of spring flow during storm events, because the simpler SCPE models adequately simulate groundwater-level and average fluid-mass-balance conditions. More recent ongoing work by Xu and others (2015) utilized a research version of CFP that includes transport between the conduit network and the SCPE (Reimann et al., 2013). It applied to a hybrid model of the Woodville Karst Plain for better simulation of long-term nitrate transport using the same long stress periods of Davis and others (2010) and improved on the HM developed by Gallegos and others (2013).

It would be useful to conduct additional tests of model approaches by comparing residence times of rapid (for example, storm event) flow components with slow flow matrix components, as derived from the use of geochemical mixing models with hydrograph separation techniques given time series data of chemistry and flow information. Unfortunately, no such time series data are currently (2022) available for the major springs within the Wakulla Springs-Leon Sinks area.