Modeling Monthly Nitrate Concentration in a Karst Spring with and without Discrete Conduit Flow

: Understanding the coupled continuum pipe ‐ flow framework for modeling contaminant transport in karst systems is critical for protecting water resources therein. This study simulated point and non ‐ point source nitrate concentration in a karst spring and investigated the results gen ‐ erated from the flow and transport model with and without discrete conduit flow. CFPv2 and CMT3D models were integrated to address the changes in nitrate concentration at a monthly scale in a karst spring, and the results were compared with that from an equivalent porous media (EPM) model with high ‐ hydraulic conductivity (K) zones set in the main karstified area to represent con ‐ duits. The results show that the CFPv2+CMT3D model is able to describe well the recession of nitrate concentration in spring discharge, and the relatively larger deviation (slower nitrate recession) from the observed trend for the EPM model is probably a result of the limitation of utilizing high hydrau ‐ lic conductivity cells to represent conduit. Moreover, simulated hydraulic heads in poorly karstified areas from the two models both show slight differences from the observations (the head RMSE val ‐ ues of calibration/validation for CFPv2 and MODFLOW models are 0.16 m/0.25 m and 0.26 m/0.17 m, respectively), indicating the inclusion of conduits may not affect the simulation considerably, and the lower the proportion of karstic area, the slight effects brought from the inclusion of conduits in the model. For highly karstified areas, the CFPv2+CMT3D model could provide more accurate results (head RMSE of calibration/validation for CFPv2 and MODFLOW are 0.22 m/0.06 and 0.52 m/0.47 m, respectively), showing the coupled continuum pipe ‐ flow framework may be more ap ‐ propriate for applying to highly and maturely karstified areas where the variations in the behavior of flow and contaminant transport are more affected by turbulent flow regime.


Introduction
Nitrate contamination in aquifers has been accelerated owing to growth in population [1] and development in urbanization [2], as well as increased requirements in the production of food [3], wastewater production, and chemical fertilizers [4].In particular, [5] studies based on a detailed global database and an elaborate analysis estimated that the number of karst water consumers in 2016 was 678 million or 9.2% of the world's population.This recent estimation is probably more realistic [6].Groundwater in karst aquifers tends to be more vulnerable to nitrate contamination due to the secondary porosity of interconnected caves, conduits, and underground drainage channels [7,8].Due to the high contaminant vulnerability, groundwater quality in karst systems is a principal management concern for water resources managers and engineers.Elevated concentration of nitrate in drinking water can cause serious health conditions (e.g., methemoglobinemia), and increased nitrogen levels may negatively affect coastal water by eutrophication which impedes the growth of other aquatic vegetation [9,10].As such, a number of techniques have been applied for predicting nitrate concentration in karst systems to help policy makers set and execute plans to control the contamination since the last century.For example, using modeling approaches, prediction for the trend of nitrate concentration could significantly help water resources managers and engineers make solid decisions on conducting groundwater quality projects and have better strategies to protect groundwater resources.
A variety of modeling approaches has been developed to investigate the groundwater flow and solute transport within the dual nature of karst systems to facilitate efficient decision making on restoration, water resources, and land use policies.Since the 1980s, numerical models have been developed and widely applied to study groundwater flow and solute transport in karst aquifers.In the past few years, the US Environmental Protection Agency (EPA) storm water management model (SWMM) flow simulator, initially developed to simulate flow in urban drainage systems, has been used as a conduit flow model for karst aquifers [11][12][13].Although the SWMM model has been employed to simulate hydraulic dynamics in previous studies, it cannot calculate solute transport in karst systems.MODFLOW-Conduit Flow Process (CFP) was developed by USGS to simulate turbulent flow in a discrete network of pipes coupled with the traditional groundwater flow equation [14], and model performance has been evaluated by [15] for a spring in Florida, in which the model performance (match between observed and simulated discharge) improved by 12-40% compared to the results from MODFLOW.Based on MOD-FLOW-CFP, a research version of CFP, named MODFLOW conduit flow process version 2 (CFPv2), has been developed by adding conduit association drainable storage (CADS, the concept of CADS considers storage volumes, where karst water is not part of the active flow system but hydraulically connected to conduits, e.g., karstic voids and large fractures.)in MODFLOW-CFP mode [16][17][18].Compared to the original version of MOD-FLOW-CFP, the CFPv2 model expands the ability to generate flow files in both matrix and conduits, which are important inputs to model solute transport.The Modular 3-Dimensional Multi-Species Transport (MT3DMS) model, which is used for the simulation of advection, dispersion, and chemical reaction of contaminants in groundwater systems, has been commonly used to solve solute transport for field-scale problems [19].The Conduit Modular 3-Dimensional Transport (CMT3D) model is developed based on the MT3DMS source code, which adds the ability to simulate solute transport in conduits to MT3DMS so that it can be integrated with the CFP or CFPv2 flow model to simulate flow and solute transport in matrix and conduits of karst aquifers [20,21].
CFPv2 and CMT3D have been integrated to successfully simulate point source nitrate transport in previous studies [22,23].However, nonpoint source nitrate concentrations in spring discharge and aquifer at a monthly scale have not been studied by the CFPv2+CMT3D model.Further, although MODFLOW has been employed to make preliminary predictions of hydrologic behaviors of a maturely karstified aquifer [24], whether the equivalent porous medium (EPM) model (MODFLOW and MT3DMS) is able to simulate groundwater flow and nitrate transport in moderately karstified aquifers is still unclear.The objective of this study is to investigate the differences between discrete-conduit and equivalent porous media models in simulating nitrate concentrations in a karst spring fed by both point and nonpoint nitrate sources.CFPv2 and CMT3D models were integrated to address the change in nitrate concentration at a monthly scale in a karst spring.The results were compared with those from the EPM model, which was developed using MODFLOW 2005 and MT3DMS.In this study, we seek to answer the following questions: (1) What are the differences between the simulation for a moderately karstified aquifer using the model with and without discrete conduit flow?(2) At a regional scale, can we use the EPM model to simulate nonpoint source nitrate-N transport for a moderately karstified aquifer?

Study Area
The Silver springshed in Florida, USA, was delineated using the potentiometric surface of the upper Floridan aquifer (UFA).The basin may undergo seasonal changes because of the temporal variation of the potentiometric head; therefore, the Silver springshed was delineated from a 1000-year capture zone and resulted in an area of approximately 2312 km 2 , as shown in Figure 1.The springshed boundary is obtained from St. Johns River Water Management District website.The springshed is spread across four north-central counties: Marion, Alachua, Putnam, and Sumter (Figure 1).Silver Springs, located in the central springshed (Figure 1), is the largest spring group in Florida and the main outlet for groundwater in the UFA, with an average discharge of 7.9 m 3 /s over the period from 2010 to 2016.The climate of the area is humid subtropical, with hot, rainy summers and cool, dry winters (Phelps 2004).The annual average temperature is 21.8 °C, and the average annual precipitation is 1290 mm (U.S. Climate data, https://www.usclimatedata.com/climate/ocala/florida/united-states/usfl0355(accessed on 20 July 2021)).The major aquifer in the study area is the Floridan aquifer, which is divided as the UFA, a middle semi-confining unit, and the lower Floridan aquifer (LFA).The Floridan aquifer is underneath the surficial aquifer and is an intermediate confining unit (Figure 2).The aquifer of particular concern is the UFA, which is mainly composed of Ocala Limestone, and the thickness is about 90 m in the study region [25,26] (Figure 2).A majority (86%) of spring discharge is from the upper 30 m of the UFA [27,28].The UFA is confined in most parts of the springshed except in some areas where the Ocala limestone crops out [28].The UFA in the western part of the springshed is near the land surface, and numerous sinkholes (more green dots appeared in the western parts of the springshed shown in Figure 1c) formed there, leading the water easily to recharge into the UFA.Therefore, the recharge rates are relatively higher, especially in the central-western portion of the springshed, where there are numerous closed sinkhole depressions [28] (Figure 1c).Subsurface karst systems, i.e., conduits, result from the dissolution of carbonate rocks, particularly by recharge through sinkhole depressions making the aquifer characterized by dual-porosity and complicating the groundwater flow and contaminant transport [28].Karst geologic features are distributed mainly around the spring outlet, where the hydraulic conductivity and groundwater flow velocity are relatively higher than those in the rock matrix [27].The horizontal and vertical positions of explored sections of water-filled conduits and caves are visible around the main spring vents, and the extent of the karst conduit system is believed to extend a significant distance from the spring vents [29].The land use and land cover (LULC) in the springshed is comprised of agriculture, residential, forest, barren land (i.e., sparse vegetation and fallow ground), rangeland, wetlands and streams, and lakes (data obtained from https://data-floridaswater.opendata.arcgis.com/(accessed on 12 February 2020)).Surface drainage in the springshed is limited (Figure 1), and the surface watershed contributes little surface runoff [27,28].Most of the rainfall either drains directly into closed depressions or seepages into the unconfined limestone of the UFA [28].As such, the transport of nitrate from land surface to spring discharge is mainly through groundwater recharge.The subsurface drainage system continues to evolve today, as evidenced in part by the frequent occurrence of sinkholes [27].There are numerous closed sinkhole depressions with permeable bottoms in the springshed (Figure 1), which are varied in dimensions and most prevalent in the mature karst terrane west of the Ocklawaha River, facilitating the recharge activity and, correspondingly, nitrate transport to the subsurface.The holes result from the collapse and flow of the sandy and clayey cover into limestone solution cavities.If a sinkhole suddenly develops near a point of groundwater discharge, such as a well or spring, temporary turbidity in the discharging water may result [27].

Nitrate Contamination to Silver Springs
Nitrate concentration in Silver Springs has increased considerably since the 1950s, and land-use practices are generally thought to affect the change in nitrate concentration in groundwater and Silver Springs [28][29][30].Nitrate in the springs could be from both nonpoint source (e.g., atmospheric deposition of nitrogen from the forested area, implementation of fertilizers in agricultural land, and residential turfgrass and golf courses) and point source pollution (e.g., wastewater leakage from onsite sewage tanks and municipal waste treatment systems as well as the nitrate-contained storm runoff with accumulated nutrient loading from impervious roads, which are mostly distributed in Marion County, especially in the western part of Silver Springs) (Figure 1) [28].Nitrogen derived from atmospheric deposition fluctuates around the mean but has remained constant since the 1950s.Fertilizers started to be used in the 1950s, and nitrogen has experienced a slight reduction since around the 1970s after a sharp increase in use.At the same time, the trend in nitrogen contamination from wastewater leakage has been upward since the 1950s due to the continuous increase in population in Marion and other counties in the springshed [28,30].In the aquifer, [28] reported that the agricultural land-use areas had the highest median nitrate concentration (1.70 mg/L) tested in a total of 30 wells, followed by residential and transportation areas (1.5 mg/L) and forested lands (0.09 mg/L).

Datasets
The observed data employed in this study include hydrologic data and nitrate concentrations in the discharge of Silver Springs and monitoring wells.The observed hydraulic heads (January 2010-December 2016, recorded at daily intervals) in the UFA in thirteen monitoring wells and nitrate concentration (December 2014-December 2016, recorded at daily intervals) in spring discharge of the main vent were both obtained from St. Johns River Water Management District (SJRWMD).The discharge record at the main vent of Silver Springs is maintained by the U.S. Geological Survey (USGS).Spring discharge is recorded at daily intervals, and the data over the periods January 2010-December 2016.The nitrate concentration data from March 2019 to March 2020 were obtained from a field experiment, which was employed to investigate the change in nitrate concentration in spring discharge.
In addition to the online datasets, a monitoring well was drilled at the UFA upgradient the springs (1.5 miles away from the springs) in March 2019 to provide observed nitrate concentrations in the aquifer.Monthly-frequency groundwater samples were collected in both the well and Silver Springs using submersible pumps from Mar 2019 to Mar 2020 and measured by Engineering, Research & Design, Inc, Orlando, FL, USA.

Model Development
Figure 3 shows the flow chart of model development for the nitrate concentration in Silver Springs.The flow model was developed using the CFPv2 model based on an equivalent porous media (EPM) model [16].Simulation of nitrate concentration in Silver Springs was accomplished by integrating the CFPv2 model with the CMT3D model [20].The details on the modeling are provided in the following sections.

The EPM Model
Model development is shown as a flowchart in Figure 3.The first step is to develop a three-dimensional groundwater flow EPM (equivalent porous media) model, which was created using MODFLOW.The model domain is discretized into 496 rows and 236 columns with a cell size of 190.5 m by 190.5 m (Figure 4a).As shown in Figure 4b, layers 1 and 2 represent the Oligocene Formation (surficial aquifer) and Miocene Formation (intermediate confining unit), respectively, according to the hydrogeologic units of the Floridan aquifer described in [27].Layer 3 represents the Avon Park Formation which presents the Ocala Limestone (i.e., the UFA).The base of the model (Layer 4) is the semi-confining layer underlying the UFA.The top of model layer 1 is at the bottom of the surface sand cover; this layer is also at the top of the intermediate confining unit.The top of model layer 2 coincides with the bottom of model layer 1 and the bottom of model layer 2 is the top of the UFA.The determination of layer thickness was based on the North Florida Southeast Georgia (NFSEG) model, which is also a MODFLOW implementation and was developed by the St. Johns River Water Management District (SJRWMD) [31].The Silver springshed is included in this model boundary.According to the data used in the NFSEG model, the thicknesses of layers 1 and 2 in the EPM model are both set uniformly to 9 m, but layer 3 can be variable, exceeding 300 m in the west but thinning to about 60 m in the central areas.According to boring logs conducted in the field near the springs, the material within the depth 0 to 3 m is mainly clayed sand.Further, [27] reported that the surficial aquifer, which is composed of marine sand and clay, ranges from 0-30 m in thickness, and the intermediate confining unit where clay and hard, dense sand were the dominant material ranges from 0-30 m.Therefore, the thickness values set in the model are within the ranges of values in both the documentation and observations.For the EPM model, the springshed (Figure 1) is set as a no-flux boundary because the springshed is delineated using a potentiometric head, and it can be considered as a subsurface drainage basin; thus, no lateral flow occurred at the boundary.The boundary conditions of the NFSEG model are a result of extensive calibration done using multiple observed groundwater wells scattered over Florida and Georgia.As such, aquifer properties and some of the boundary conditions (e.g., rivers and lakes) for layers in the domain of the EPM model remain the same as those used in the NFSEG model.Specifically, a head-dependent flux boundary was used to simulate rivers, lakes, as well as streams; the boundary condition for these water bodies was determined according to the NFSEG model.Only one discharge point in the springshed, Silver Springs, is assumed for the model and treated as a drain boundary.
Hydraulic conductivities from the NFSEG model were used in the EPM model.As described, [31] used groundwater heads at monitoring wells across the domain to calibrate the hydraulic conductivities in the NFSEG model.In the EPM model, calibrated horizontal hydraulic conductivities in layer 1 range from 3 × 10 −5 m/s to 5 × 10 −5 m/s at the central and western parts of the study region.In other areas, conductivities range from 4 × 10 −6 m/s to 3 × 10 −5 m/s.The horizontal hydraulic conductivity distribution in layer 2 is similar to but much higher than layer 1 at the central and western parts of the study area.Horizontal hydraulic conductivities range from 3 × 10 −4 m/s to 0.06 m/s caused by the cropping out of Ocala limestone in some parts around Silver Springs.Hydraulic conductivities in layer 3 are much higher than that in layer 2 (e.g., 0.02 to 0.3 m/s in the central and western parts of the study area) due to fractures and conduits, but the distribution is similar to layer 2. Vertical hydraulic conductivities in layers 1 and 2 are much lower than horizontal conductivities, in particular in the north and south of the study region, because the aquifer consists of a low permeable rock matrix with clay or other horizontal low permeability layers.For reasons discussed earlier, conduit networks are assumed to be predominantly in the deeper part of the aquifer and are simulated in layer 3.According to [31], high hydraulic conductivity cells were used for representing karst features in the UFA in the EPM model (Figure 5a).

The CFPv2 Flow Modeling
CFPv2 simulates flows in both rock matrix and conduits.Darcy's law is applicable in the rock matrix where the flow is laminar, and the governing equation is [33]: where  ,  , and  (LT −1 ) are the hydraulic conductivity along the -, and  directions, respectively; ℎ(L) is the hydraulic head;  is the source and/or sink term (T −1 );  (L −1 ) is the specific storage of the porous material; and (T) is time.However, Darcy's law is only appropriate for use in laminar flow condition.The Darcy-Weisbach equation for a pipe is applicable for non-laminar groundwater flow in a conduit [14], which can be expressed as: where ∆ℎ and ℎ are head loss along pipe length ∆(L),  is the friction factor,  is the pipe diameter (L),  is the mean velocity (LT −1 ), and  is the gravitational acceleration constant (LT −2 ).Caves and conduits are defined as a pipe network consisting of discrete cylindrical tubes and several nodes in the model.The nodes located in each MODFLOW cell compute the exchange between tubes and porous matrix.The linear equation used by CFPv2 to calculate flow exchange between the porous matrix and a conduit node is given by where  is the volumetric exchange flow rate (L 3 T −1 ),  , , is the pipe conductance at cell , ,  (L 2 T −1 ), ℎ is the head at conduit node , and ℎ , , is the head at cell , , .

The CMT3D Modeling
Similar to CFPv2 adds the ability to simulate flow in conduits to MODFLOW, CMT3D equipped with the ability to compute contaminant transport in conduits to MT3DMS [19].The general three-dimensional transport equation used in MT3DMS is applied to simulate solute transport in matrix of porous medium.In this study, only one species solute nitrate-N is considered, where  (−) is the porosity of porous medium,  (ML −3 ) is the nitrate concentration,  (T) is time,  (L) is the distance along the i direction,  (L 2 T −1 ) is the hydrodynamic dispersion coefficient in porous medium,  (LT −1 ) is the linear pore water velocity,  (T −1 ) is the volumetric flow rate per unit volume of aquifer representing sources and sinks,  (ML −3 ) is the nitrate concentration of the source and sink flux, and  (ML −3 T −1 ) is the chemical reaction term which is not considered in this study.
MODFLOW and MT3DMS cannot simulate flow regimes, and solute transport in karst features well using high hydraulic conductivity cells and is limited to laminar conditions.Therefore, the CFPv2 and CMT3D model was used to simulate flow and contaminant transport in the Silver springshed.The mechanisms for the two models are described in the following sections.
Other than MT3DMS, which can only simulate contaminant transport in a porous medium, CMT3D is able to model contaminant transport in discrete conduits employing flow output from CFPv2 [16].CFPv2 produces the link file to MT3DMS for matrix transport and the link file to CMT3D for conduit transport.Subsequently, CMT3D uses these files to compute transport in matrix using MT3DMS, then calculating transport in conduits using an advective transport module [34]: where  (ML −3 ) is the nitrate concentration in conduit tube,  (L 2 T −1 ) is the dispersion coefficient.

Conversion of EPM to CFPv2
The CFPv2 model was developed based on the EPM model, and the model boundary conditions described above were retained in the CFPv2 model.The steady-state CFPv2 model was created by incorporating conduit networks which were simulated by conduit node and discrete cylindrical pipes [14], into the UFA (the third layer of the model), where karstic features were mostly distributed [27,28].Ref. [27] delineated a fracture and possible faults map in central Florida.However, it is still a challenge to identify conduit structures in this area [35], in particular in Silver springshed.Conduit structures in this study were determined based on sinkhole locations.Since groundwater enters the system both through surface infiltration and direct injection of runoff into karstic conduits via sinkholes, locations of the sinkholes provide an indication of where the conduits are located [8,27].Therefore, the locations of conduits were assigned according to the distribution of sinkholes, i.e., the model uses known sinkhole locations to characterize the karst arrangement, as seen in Figure 5b.Specifically, the cells where a high density of sinkholes around were replaced with CFPv2 pipes.Accordingly, there are seven conduits developed in the CFPv2 model (Figure 5b), which are comprised of 687 discrete cylindrical conduit pipes.As described in Section 2.4.2, the cells which were given high-K values in the EPM model were represented as conduits.Here some of these cells were replaced with discrete conduits; and because hydraulic conductivity in the matrix around the conduits could be low (i.e., lower than their K values in the EPM model), the hydraulic conductivity values in these cells were in succession of lowering to the average hydraulic conductivity values of the surrounding matrix [22,32].A tracer study conducted by [29] reported that there were no mature conduits in the upper portion of the springshed; therefore, conduits were assigned mostly in the central-west portion of the domain (Figure 5b).The uncertainty of conduit structure was not considered in this study.Ref. [36] established an approach to evaluate spatiotemporal dynamics of controlling parameters in distributed environmental models.Future research may use this approach to test the sensitivity of conduit structure to modeling results.The boundary conditions of the EPM model described in Section 2.4.2 were retained in the CFPv2 model.Recharge values were assigned variable spatially across the springshed based on the NFSEG model [31].These recharge values were used as the long-term average recharge values from January 2010-December 2016 by specifying the recharge factor as one in the MODFLOW-CFP program for the steady-state simulation.The steady-state model used a single, one-day time step.
In order to investigate and project the annual nitrate concentration in Silver Springs, a transient flow simulation was conducted from January 2010-December 2016.Monthly recharge values (recharge factors in the MODFLOW-CFP program) were quantified by calibration using the observed hydraulic heads and spring discharge in the main vent in the transient simulation from January 2010-December 2016.Monthly stress periods with a daily time step were used.Effective porosity governs the velocity of groundwater movement in transient simulations.Fine to coarse-grained sands are the dominant material in the surficial aquifer; therefore, the simulated effective porosity for layer 1 is set as 0.25.Model layer 2 is a semi-confining layer where clay is the main soil type combined with sands and some hard limestone so that the simulated effective porosity for layer 2 is set as 0.05.An effective porosity of 0.3 was used for layer 3 because the aquifer was loose and fractures were highly distributed.These porosity values are also within the ranges of porosity documented in [28] (0.15-0.4).

Development of CMT3D Model
When groundwater flow is simulated by the CFPv2 model, the transport model in this study, the CMT3D model, will use the hydraulic head and flow data in both the rock matrix and conduits generated by the transient CFPv2 model to simulate nitrate transport in matrix and conduits.Figure 3 shows the basic integration mechanism of CFPv2 and CMT3D.CFPv2 generates the link-MT3DMS file for matrix transport (LMT6) and the link-CMT3D file for conduit transport (LCT6).CMT3D uses these files to sequentially calculate transport through the matrix using MT3DMS and then calculate transport through conduits using an advective transport model [37].The flux between the matrix and conduits is calculated using Equation (3).
Nitrate reduction is negligible in the UFA and other aquifers due to low ammonium and dissolved organic carbon (DOC) concentrations, and significant nitrate reduction is unexpected when there is low carbon availability [38][39][40][41][42][43].Moreover, the flow velocity is high in conduit, especially around the Silver Springs, where karst features are maturely developed, so the residence time for the solute could be short.Moreover, conduit transport dispersion and chemical reaction were not considered in the simulation of contaminant transport in the Woodvill Karst Plain, which is also located underneath the UFA, due to the limited residence time [22].Therefore, nitrate is set as the non-reactive solute in this study, and conduit transport dispersion and chemical reaction are not considered.
In the transport model, three parameters are required to be solved to simulate nitrate concentration in Silver Springs: dispersion coefficient in the rock matrix and conduits and nitrate concentration in groundwater recharge (i.e., input nitrate concentration) [32].A tracer test was conducted to estimate the dispersivity.The tracer dye (rhodamine WT) was injected into a monitoring well (new drilled well marked in Figure 1c) installed in the UFA near Silver Springs.The tracer concentration was sampled using Sigma 900 MAX portable sampler at a daily frequency and tested from May 2019 to March 2020.The two-region non-equilibrium transport model in the program CXTFIT [44] was applied to analyze the tracer test dataset [45].Nitrate concentration in groundwater recharge, which is assigned as a constant head boundary in the transport model, is described in the following section.The study period of the transport model was the same as that of the transient flow model, which was from January 2010-December 2016 and monthly stress periods with a daily time step were used.

Nitrate Concentration in Groundwater Recharge
Nitrate-N concentration in groundwater recharge from nonpoint (fertilizer in agricultural and deposition in forested lands) and point source (septic tanks in residential and storm runoff in transportation areas) are both modeled via land use and land cover.In this study, the nitrate concentration in groundwater recharge ( , [M/L 3 ]) were assumed to be linearly proportional to population density ( [pop./L 2 ]): The slope  [M/L 3 ] and the background concentration  [M/L 3 ] were estimated during the calibration of the transport model.
The details for the estimation of nitrate-N concentration in groundwater recharge were described in [32].The nitrate-N concentration in groundwater recharge was set as 4 mg/L for agricultural land and 0.12 mg/L for forested land [28].As population increases, the nitrate-N concentration leaked from septic tanks may also increase within a given residential area.Further, more impervious areas (i.e., transportation corridors) may be developed as a result of population growth, and more stormwater runoff from these areas may be generated.Therefore, for residential and transportation lands, the nitrate concentration in groundwater recharge was assumed to be linearly proportional to population density [32].Parameters of the linear function (slope and intercept) were estimated during the calibration of the transport model.

Model Calibration and Validation
The trial-and-error method was used for calibration in this study.Model calibration begins with hydraulic calibration of the steady-state model, recharge calibration of the transient model, then progresses to calibration of the nitrate transport model.Hydraulic conductivity, conduit diameter (), tortuosity, and conductance ( , , ) were calibrated in the steady-state CFPv2 model.Calibration and validation targets include observed water levels in 13 monitoring wells (Figure 1) and observed discharge of the Silver Springs.Periods for calibration and validation are from January 2010 to December 2013 and from January 2014 to December 2016, respectively.Calibration of conduit parameters was conducted by comparing the mean observed and calculated water levels in 13 monitoring wells and discharge of Silver Springs over the period from 2010-2016.Particularly, the hydraulic conductivity values around the conduits were lowered to the average hydraulic conductivity values for the surrounding matrix after calibration.Calibration of recharge is required for the transient model.The recharges of the entire model across the model boundary come from several sources: (1) Net precipitation (rainfall minus evapotranspiration), (2) Creek flow into sinkholes, (3) Irrigation at farm lands, and (4) Discharge from onsite sewage disposal system.The calibration applied adjustment factors to recharge estimates for use in the model.Finally, nitrate concentration in groundwater recharge for each source was calibrated based on the observed annual average nitrate concentration in spring discharge during December 2014-December 2015, and the nitrate transport model was validated during January 2016-December 2016.

Steady-State Models
The steady-state models were calibrated based on the observed groundwater levels and spring discharge.Calibrated values for the conduit diameter () and pipe conductance ( , , ) varied throughout the model domain.As shown in Table 1, calibrated values for conduit diameter range from 3-6 m, with higher values representing the main west to east running conduits that feed Silver Springs.Calibrated values for pipe conductance varied from 9.3 10 to 9.3 10 m 2 /day, with higher values indicating increased flux between conduits and the surrounding matrix.Tortuosity was calibrated with a final value of 1.5 for all the conduits.Although these parameters cannot be measured directly, the results are within previously reported ranges for karst aquifers.For example, conduit diameter in the UFA was documented as ranging from less than 1 m to 15 m [22], and pipe conductance values were allowed to vary between 3.0 10 m 2 /day and 1.8 m 2 /day [46].Ref. [47] characterized channel networks in carbonate aquifers and found that tortuosity was within the range of 1.103 to 2.369.Similar results were documented in [48], in which the tortuosity values were reported ranging 1 to 3.9.As shown in Figure 6a, the observed values versus simulated steady-state heads from CFPv2 model accumulated on the 1:1 line, indicating that the modeling results match the observed counterparts.The calibration statistics for the steady-state models are shown in Table 2.The average root mean square error (RMSE) and relative error between simulated and observed potentiometric heads across monitoring wells (Figure 1) are 0.28 m and 2%, respectively, both of which are lower than that of the EPM model (0.36 m and 3%).These statistics indicate that the steady-state CFPv2 model is able to generate more accurate results in comparison with the EPM model.The calibrated recharge adjustment factor for the transient hydrologic models considers both wet and dry months.Based on the recharge estimated in the NFSEG model, the adjustment factor corresponding to the recharge in the steady-state simulation (the first stress period) was assigned as 1 in the transient model.The adjustment factor ranging from 1.5-3 was applied in the wet months with higher recharge, while values from 0.2-0.9 were used for the months with lower recharge, and values from 0.9-1.5 were employed in the normal months.
The scatter plot of observed versus simulated heads for the transient CFPv2 simulation shows a good agreement with no apparent over-or under-estimation bias for both calibration and validation (Figure 6b).The RMSE of simulated heads for calibration and validation are 0.33 m and 0.29 m, respectively, and the relative error is 1.8% (1.6%) for the calibration (validation) period.Both statistics are lower than that of the EPM model (Table 2).Figure 7 shows the change in observed and simulated hydraulic heads throughout the simulation periods in three wells located in different regions of the springshed (Figure 1).Well 10201167 is located in the upper part of the springshed, where karst features are sparsely distributed; thus, the karst aquifer is not well developed.The similar settings of hydraulic conductivity values in this region cause no expressive differences between the simulation results from the two models (Figure 6a), in which the head RMSE values of calibration/validation for CFPv2 and MODFLOW models are 0.16 m/0.25 m and 0.26 m/0.17 m, respectively.For well 70232656, which is an approach to the central part of the springshed where karstic features are mainly distributed [28], the simulation differences between the models are enlarged and the simulated heads from the CFPv2 model are much closer to the observed ones (Figure 7b; head RMSE of calibration/validation for CFPv2 and MODFLOW are 0.22 m/0.06 and 0.52 m/0.47 m, respectively).The results indicate that although there are also a few karstic features found in this region (Figure 1), thus no evident differences between the settings of models, the inclusion of the conduits in the central part of the springshed can influence the simulated head values in the area nearby.When the well is located adjacent to the conduit (well 30272814; Figure 1), the simulated hydraulic heads from the CFPv2 model are much closer to the observed counterparts (Figure 7c; head RMSE of calibration/validation are 0.03 m/0.02 m and 0.12 m/0.11 m, respectively, for CFPv2 and MODFLOW models).Particularly, probably due to the high groundwater flow velocity in the conduit decreasing the recession of the potentiometric head, the simulation from the EPM model tends to have a longer recession of head values compared to that from the CFPv2 model, although it is capable of capturing the main trend of the heads.Further, Figure 8 shows the simulated discharge of Silver Springs from the two models are both able to capture the spring hydrograph over the study period; however, the CFPv2 model creates a much closer agreement to the observed hydrograph.Similar to the simulated potentiometric heads in well 30272814 (Figure 6c), the inclusion of the conduits in the CFPv2 model has the groundwater flowing faster, thus creating steeper spring recession curves than those from the EPM model.A further indication of successful calibration and validation is the high Nash-Sutcliffe efficiency (NSE) values for the spring discharge, which are 0.85 (calibration) and 0.84 (validation), respectively (Table 1).The performance of the transient CFPv2 hydrologic model is better than that of the EPM model, as indicated by the statistics in Table 2 and Figure 8, showing that the inclusion of conduits instead of high hydraulic conductivity cells can better compute non-laminar flow in the karstic medium.

Nitrate Transport
Longitudinal dispersivity in the rock matrix yielded from the tracer data is 1.8 m.The ratio of longitudinal to transverse dispersivity ranged from 0.05 to 0.1, with the upper bound set based on reported maximum expected values [19].The ratio of longitudinal to transverse dispersivity is set to 0.1.Both values are on the same magnitudes as previously reported ones in the UFA [49,50].
Model performance at springs is illustrated in Figure 9.In some periods, the observed nitrate concentration in spring discharge generally follows the trend of observed spring discharge, where concentration reached high values at peaked spring discharge (Feburary 2015) and decreased gradually as the discharge reduced (Feburary 2015-July 2015, Figure 8b).The observed nitrate concentration does not follow the trend of spring discharge in other periods (April 2015-July 2015, August 2016-December 2016) because the recession of nitrate is slower than that of spring discharge.The relation between hydrologic conditions and concentration is similar to that reported in [51], in which the observations of nitrate time series over multiple seasons were investigated, and it shows the higher concentrations coupled with higher discharge values contributed to higher nitrate fluxes in the higher hydrologic conditions.Changes in simulated nitrate concentration in spring discharge show weak fluctuations from December 2014 to December 2016, and both of the simulations can capture the main trend of the observed concentration, except for the period from December 2015 to April 2016.In this study, nitrate in the springs is derived from both point and nonpoint sources, which are mostly distributed in non-or slightly karstified areas (agricultural and forested lands) that account for a relatively larger part of the springshed [32].The slow recession of nitrate in spring discharge is probably due to the slow solute transport from the nonpoint source to the springs via a low hydraulic conductivity medium.The factors that resulted in the change in nitrate concentration not following the spring discharge may also come from other sources besides the change in LULC and the growth of the population.For example, ref. [52] stated that if nitrate concentrations in the soil and epikarst were high prior to storm events, fast increases in nitrate concentrations can be observed as storm responses; if nitrate concentrations in the soil and epikarst were low prior to storm events, heavy rainwater may dilute nitrate concentration in spring water.Moreover, heavier fertilizers might be applied to non-agricultural lands (e.g., recreational land and home lawns) [53], which may be a considerable contribution of nitrate accumulations in soils.Correspondingly, the behavior of slow recessions of nitrate concentration may be due to a small amount of nitrogen in the soil prior to normal or even larger storm events that cause the dilution effect of nitrate concentration increase.For the CFPv2+CMT3D model, the RMSE for the calibration (validation) period is 0.02 mg/L (0.11 mg/L), and the relative error for the calibration (validation) period is 1.5% (8.3%).The model performance improved compared to the EPM model, indicated by larger statistics where RMSE for the calibration (validation) is 0.04 mg/L (0.12 mg/L), and the relative error for the calibration (validation) is 3% (8.5%).The data used for calibration and validation for hydraulics and transport are large enough.Hence, the results of the decreased simulation error from the EPM model to the CFPv2+CMT3D model are solid.The model performance was also tested by Saller et al. (2013), in which the match to observation well hydraulic heads was substantially improved with the addition of conduits in the Madison aquifer, South Dakota, USA.The reason why the simulation has a relatively small improvement is that most part of the aquifer within the springshed has not been well karstified.Since the reaction module has not been included in the CMT3D model by far, more data to test the performance of the transport model and whether decreased error from use of the EMP model to the CFPv2+CMT3D model could offset model complexity (e.g., inclusion of reaction module) will be considered as future research.The relatively larger deviation (slower nitrate recession) from the observed trend for the EPM model probably results from the limitation of utilizing high hydraulic conductivity cells to represent conduit; thus, the turbulent flow may not be counted.Further, both of the models overestimate the nitrate concentration from November 2015 to April 2016.The reduction in nitrate over this period is probably due to the spring restoration projects that had been undertaken to reduce nitrogen loading from human activities to protect Silver Springs, especially from 2015 to 2016 [54].
There are several broad observations that can be gleaned from the model performances in 2019 (Figure 9b).First, slight differences exist between the two simulated nitrate concentrations in spring discharge (RMSE are 0.08 mg/L (CFPv2+CMT3D) and 0.09 mg/L (EPM), and relative errors are 6% (CFPv2+CMT3D) and 7% (EPM)) except in the end of period (Figure 9a; October 2019-December 2019).The divergence of nitrate concentration and spring discharge during the end of the period (October 2019-December 2019) could also be due to the slow recession of nitrate concentration.The trend of nitrate concentration follows that of the spring discharge until October 2016.After then, the spring discharge decreased while nitrate concentration still climbed slowly.Since a large part of the aquifer within the springshed has not been maturely karstified (upper and lower parts), the velocity of the nitrate source transports to the springs is still low.Although the spring discharge underwent a decreasing trend by small storm events that occurred in the region where the springs were located, the nitrate source could still transport to the springs in the same period from other parts of the springshed.This causes the slow response of nitrate concentration to spring discharge.Similar to the previous analysis, the simulated trend follows the spring discharge and results from the EPM model tends to have a slow recession relative to that from CMT3D.
A second observation is that there is a distinct pattern of simulated nitrate trend compared to the observations in 2019 (Figure 9b).During this period, variations in nitrate did not respond to the hydrologic conditions.Nitrate increased from April to June during the period when the discharge reduced, whereas discharge gradually decreased from September, corresponding to an increase in nitrate.The factors that resulted in the change in nitrate concentration have been stated previously [52,53].The behavior of peaked and the increasing trend of nitrate concentration may be due to the accumulated large amount of nitrogen in the soil prior to normal or even small storm events that increased the nitrate flux, and vice versa for the low nitrate data in this period.Ref. [52] studied the various nitrate responses in karst aquifers to storm events and nitrate patterns at the springs, which elucidates the relation of nitrate responses at karst springs with hydrologic conditions and nitrate availability.Future research may consider the relation between nitrate loading in soil and aquifer, the mechanisms of the response of storm events to transport of nitrate from soil to karst aquifers, and the effects of storm events on nitrate concentration in spring discharge.If these factors are included in the model, more solid simulation results would be obtained.
Finally, the simulated nitrate in groundwater (new drilled well shown in Figure 1c) reveals that the results from the CFPv2+CMT3D model are lower compared to the EPM model and have a closer agreement with the observations (Figure 9c).Since found sinkholes in the springshed have been repaired, transport of nitrate from both point and nonpoint sources into the UFA could be through soil and epikarst.Even though septic tanks in the City of Ocala (point source) are set up underground, the tanks may not be put in conduits during the construction.Therefore, in the CFPv2+CMT3D model, nitrate from both point and nonpoint sources is set as no direct recharge to the conduits.Hence, nitrate in the conduit can be only from the surrounding rock matrix.When the heads in the matrix are higher than those in conduits, groundwater may flow from the matrix to the conduits [16,17].In succession, groundwater with high flow velocity in the conduits delivering nitrate to the springs may further decrease the concentration in the surrounding matrix while there is no such exchange with respect to groundwater and solute, and this is the reason why the simulated nitrate concentrations in the aquifer are higher in comparison with the coupled CFPv2 and CMT3D model.
The findings described above may be relevant to the karstified degree of an aquifer.Although well-developed conduits and caves are found around the main spring vents [29], the proportion of karstic features in the springshed is not big, and most of them are located in the central springshed [27] (Figure 1c).As such, the inclusion of conduits may not affect the simulation considerably, indicating that the lower the proportion of karstic area, the slight effects brought from the inclusion of conduits in the model.The simulated hydraulic head for well 10201167 located in non-karstified region could be an illustration for this point (Figure 7a).It is not unreasonable to consider that the variation in the behavior of nitrate pattern may be more dramatic for a higher karstified aquifer; for example, conduits are also widely distributed in other parts of the springshed, in which turbulent conditions may account for a larger proportion of flow.In this case, the inclusion of conduits by the coupled CFPv2 and CMT3D model may be more appropriate for simulating flow and transport than using high hydraulic conductivity cells to represent karstic features.

Sensitivity Anaysis
Sensitivity analysis was conducted to evaluate the influence of variations of conduit properties and solute transport parameters (porosity) on simulation.To study the impacts of variation of conduit parameters on the simulated spring discharge, each conduit parameter (diameter, pipe conductance, or roughness height) was increased/decreased from the calibrated values by 25%, 50%, and 75%, respectively, while all other parameters keep unchanged.
Figure 10a shows the sensitivity of three parameters for spring discharge.The most sensitive parameter is pipe diameter.By increasing (decreasing) the diameter by 25%, spring discharge may have a 25% growth (reduction), and increasing (decreasing) the diameter resulted in an increase (decrease) in spring discharge.The sensitivity of spring discharge on pipe conductance is weaker compared to that of diameter, and increasing (decreasing) the conductance causes growth (reduction) of spring discharge.In this study, no direct recharge to the conduits was set up in the simulation; therefore, water in the conduits is from the surrounding matrix, and increased water exchange between conduit and matrix may result in a higher volume of water flows into the conduits so that much water can be delivered to the spring.Pipe roughness height is the least sensitive of conduit parameters.Less than 2% increase (decrease) in spring discharge was caused by a 25% change in the parameter.The impact of porosity on the simulation of nitrate-N concentration in spring discharge was conducted by increasing the calibrated porosity value (0.3) in layer 3 to 0.6 or decreasing it to 0.03.The increase of porosity results in lower nitrate-N concentration in spring discharge, and decreased porosity leads to higher variation in nitrate-N concentration (Figure 10b).When the porosity decreases, more contaminants tend to be stagnant in

Change of Spring Discharge (%)
the matrix with higher capillary barriers and lower flow velocity [55].Therefore, the response of concentration to hydrologic conditions, i.e., recharge, could be more sensitive than that in a higher porosity matrix that has lower stagnant contaminant and higher transport capability.

Conclusions
In this study, an integrated CFPv2 and CMT3D model was employed to simulate point, and nonpoint source nitrate-N concentration in the discharge of Silver Springs, FL, and the simulation was able to capture the monthly trend of nitrate in spring discharge successfully.The results were compared with that from an EPM model with high-K zones set in the main karstified area to represent conduits.In the former, modeled conduit locations were based on the distribution of sinkholes.Groundwater flow input for the CMT3D model was supplied by the CFPv2 model, and input nitrate-N concentration to the CMT3D was estimated based on land use change and population growth.The flow models were evaluated using hydraulic heads in observation wells and spring discharge, while the transport model was evaluated by nitrate-N concentration in spring discharge.
Hydrologic calibration results indicate that the CFPv2 model more accurately represents the discharge in the Silver Springs and head values throughout the aquifer.Although both simulated nitrate-N concentrations from the two models are able to describe the observations, the results acquired improvement by the integrated CFPv2 and CMT3D models.Slower nitrate recessions were observed in the simulation from the EPM model, which probably resulted from the limitation of utilizing high hydraulic conductivity cells to represent conduit, and the turbulent flow may not be counted in the calculation.Further, the simulated nitrate-N concentrations in groundwater by the integrated model generated a much closer agreement with the observations in a monitoring well near the springs, indicating the inclusion of conduits in the model may improve the simulation since processes of groundwater and solute exchange between the conduits and the matrix are considered in the calculation.In addition, simulated hydraulic heads in poorly karstified areas in the study region from the two models both show slight differences with the observations, whereas the integrated model could provide more accurate results in highly karstified areas, indicating that the coupled continuum pipe-flow framework may be more appropriate for applying to highly and maturely karstified areas where the variations in the behavior of flow and contaminant transport are more affected by turbulent flow regime.
Although the general trend of nitrate in spring discharge is well simulated by the model in this study, the variations in spring discharge and nitrate (e.g., dilution) response to storm events in the model should be considered.Further research could focus on the simulation of responses of nitrate-N concentration at a daily scale in a karst spring with high and low nitrate concentration in the soil to storm events using models.Chemical reactions that take place during transport may also be considered in future research to improve the modeling accuracy.Results of such studies will allow for improved simulation of solute transport in karst springs and aquifers.

Figure 1 .
Figure 1.(a) The overall view of the U.S. and the location of the Silver springshed; (b) location of the Silver springshed in the state of Florida, USA; (c) the observation wells, sinkholes, Silver Springs, and Silver River within the springshed; well 10201167 and well 70232656 are located in the upper part of the springshed (slightly karstified areas); the new drilled well and well 30272814 are located in the central part of the springshed (highly karstified areas).

Figure 3 .
Figure 3.The procedures for modeling nitrate concentration in Silver Springs.

Figure 5 .
Figure 5. (a) Hydraulic conductivity values set in layer 3 of the EPM model [31]; (b) the delineated conduits in the UFA and the distribution of sinkhole locations in the springshed.

Figure 6 .
Figure 6.(a) Observed long-term average water levels versus the simulated hydraulic heads from the steady-state CFPv2 model; (b) observed water levels versus the simulated hydraulic heads from the transient CFPv2 model.

Figure 7 .
Figure 7.The observed and simulated potentiometric heads in (a) well 10201167 (located in the upper poorly-karstified region of the springshed), (b) well 70232656 (located approach to wellkarstified region of the springshed), and (c) well 30272814 (located in the central conduit-adjacent region of the springshed).Well numbers and locations are shown in Figure 1c.

Figure 8 .
Figure 8.The observed and simulated discharge of Silver Springs.

Figure 9 .
Figure 9. Simulated nitrate concentration in the discharge of Silver Springs for the periods (a) January 2014-December 2016 and (b) March 2019-December 2019; (c) simulated nitrate concentration in a monitoring well (new drilled well shown in Figure 1c) in City of Ocala, FL, USA.December 2014-December 2015 and January 2016-December 2016 are the calibration and validation periods for the transport model, respectively.Nitrate concentrations were both sampled at the discharge of Silver Sprigs main vent and the new drilled well (shown in Figure 1c) from March 2019-December 2019 to investigate the effects of conduit flow to model performance.

Figure 10 .
Figure 10.(a) Discharge sensitivity to pipe diameter, pipe conductance, and roughness height; (b) sensitivity of nitrate-N concentration in spring discharge to porosity.

Table 1 .
The calibrated values of diameter, conductance, roughness height, tortuosity for each conduit.

Table 2 .
Hydrologic calibration statistics for the EPM and CFPv2 model.