Modeling Downward Groundwater Leakage Rate to Evaluate the Relative Probability of Sinkhole Development at an Under-Construction Expressway and Its Vicinity

Sinkhole development has been recognized as a major geohazard, as sinkholes pose great threats to infrastructure, such as buildings, roads, bridges, and pipelines, resulting in huge financial losses to society. Previous studies indicated that the spatial density of sinkholes increases linearly with the downward groundwater leakage rate (DGLR) (inter-aquifer flow rate from an unconfined to a confined aquifer through the aquitard between them) and that the spatial variation of annual-average DGLR is a useful indicator of the relative probability of sinkhole development. In this study, a groundwater flow model using the MODFLOW computer code was developed and calibrated to simulate the spatial variation of annual-average DGLR to evaluate the relative probability of sinkhole development at an under-construction expressway and its vicinity. The results indicated that the expressway construction site has a relatively high probability of sinkhole development in the designed range of the pavement structures, and it is concluded that engineering action should be taken in advance to minimize potential sinkhole hazards.


INTRODUCTION
The major landforms known as cover-collapse and cover-suffosion sinkholes (abbreviated as "sinkholes" in the following text) have been identified as the primary geohazards in central Florida (United States) since the 1950s and have resulted in large financial losses to society, especially in populated cities, due to the catastrophic damage they cause to buildings, roads, bridges, and pipelines (Maroney et al., 2005;Lindsey et al., 2010;Brinkmann, 2013;Kuniansky et al., 2015). Sinkholes were identified statewide by field surveys and geomorphological mapping from 1950 to 2014 and are recorded in a downloadable subsidence incident report database created by the Florida Geological Survey (one of the divisions of Florida Department of Environmental Protection). In spite of drawbacks such as under-reporting of sinkholes in rural areas, the subsidence incident report database has been validated to be the most complete, accurate, and representative sinkhole inventory statewide Fleury et al., 2008), and previous studies regarding sinkhole research in Florida have been developed through the use of these data.
In central Florida mantled (buried) karst terrains, the reported sinkholes are not evenly distributed: some areas have many, while others have none (Gray, 2014). The underlying cause of sinkhole formation is limestone dissolution/removal by infiltrated weakly acidic rainwater (dissolution/removal rate extremely slow; as little as millimeters per thousand years), and the major triggering factors include heavy rainfall, rapid decline of potentiometric level due to groundwater overexploitation, and high rate of downward groundwater leakage (Beck, 1986;Ford and Williams, 2007;Brinkmann and Parise, 2009;Gutiérrez et al., 2014). Here, downward groundwater leakage rate (DGLR) refers to inter-aquifer flow rate from an overlying unconfined aquifer to an underlying confined aquifer through an aquitard (see the black arrow in Figure 1). Jammal (1982) found that the sinkholes reported in Winter Park in central Florida are more likely to form at the beginning of the wet season when DGLR rises to its annual highs. Wilson and Beck (1992) pointed out that the reported uneven distribution of sinkholes in the greater Orlando area in central Florida might be attributed to spatial variation in annualaverage DGLR. Xiao et al. (2016Xiao et al. ( , 2018 conducted statistical analyses of the relationships between the spatial distribution of reported sinkholes and the spatial variation of annualaverage DGLR in central Florida and found that sinkholes are likely to occur in those areas that have higher DGLR since sinkhole spatial density increases linearly with DGLR. Based on these findings, spatial variation of annual-average DGLR has been recognized as a useful indicator of the relative probability of sinkhole development in central Florida. The main limitation to quantifying it analytically in the central Florida sinkhole-prone areas is a lack of observed data, since local-scale groundwater observation systems have not been installed ubiquitously, and monitoring work has not been conducted routinely. Due to the recent rapid development of high-speed computers and simulation codes (Zhou and Li, 2011;Anderson et al., 2015), however, groundwater models have been successfully applied in many case studies worldwide for quantifying DGLR numerically (Sanford, 2002;Scanlon et al., 2002;Chitsazan and Movahedian, 2015;Mali and Singh, 2016;Sahoo and Jha, 2017), and these can be adopted to evaluate the relative probability of sinkhole development in central Florida.
The objective of this study is to evaluate the relative probability of sinkhole development at an under-construction expressway and its vicinity in northern Orlando, central Florida (United States), based on spatial variation of annual-average DGLR quantified by developing a groundwater model using the MODFLOW computer code, aiming at (1) providing a good reference for determination of whether the expressway construction site and its vicinity are at risk for sinkhole development, (2) providing a warning that immediate action (from engineering perspectives) should be taken in advance to minimize the hazards of sinkhole development at the designed range of the pavement structures, and (3) providing knowledge and understanding of sinkhole development in central Florida mantled karst terrains for engineers/technicians working in the field of sinkhole hazards, as well as for residents dwelling in sinkhole-prone areas.

OVERVIEW OF THE STUDY AREA
The study area, shown in Figure 1, is an under-construction expressway (see the solid white line) that connects two important state roads (SR46 and SR 429) (see the black dot-dash line) and its vicinity in northern Orlando (see the red star), central Florida (United States). The boundaries of the area are either parallel or perpendicular to the water table contours generated by the output from the north-central Florida groundwater flow model (St. Johns River Water Management District), and regional groundwater flow directions are from southwest to northeast and southeast (see the light green arrows). The land surface elevations vary from 5 to 30 m [NAVD88], with a regional average of approximately 16-17 m [NAVD88]. The study area has a subtropical and humid climate with humid/hot summers (mean maximum temperatures exceeding 30 • C) and dry/mild winters (mean minimum temperatures dropping below 10 • C) with mean annual rainfall varying from 1200 to 1300 mm (wet season from June to October and dry season from November to May) (Tibbals, 1990). The hydrostratigraphic units shown in Figure 1 consist of (from top to bottom) the (unconfined) surficial aquifer, the upper confining unit (aquitard), the (confined) Floridan aquifer, and the lower confining unit (aquiclude); detailed descriptions of their characteristics can be found in Miller (1986) and Williams and Kuniansky (2016). The surficial aquifer (thickness varying from 5 to 15 m) is primarily composed of medium to medium dense sand with moderate transmissivity. The Floridan aquifer (thickness varying from 600 to 650 m) is primarily composed of limestone/dolomite with high transmissivity (varying from 500 to 100,000 m 2 /day) since it is highly karstified, with sinkholes, sinking streams, and springs present over most of its extent. The upper confining unit (thickness varying from 20 to 25 m) is primarily composed of clay, silty clay, and sandy clay with relatively low permeability, which limits inter-aquifer flow from/to the overlying surficial aquifer to/from the underlying Floridan aquifer.

NUMERICAL MODELING
A groundwater model was developed to quantify spatial variation in annual-average DGLR to evaluate the relative probability of sinkhole development at the under-construction expressway and its vicinity, and the groundwater model developed was calibrated using the trial-and-error method by adjusting the values of horizontal/vertical hydraulic conductivities until the simulated groundwater levels and spring discharges matched the fieldmeasurements to a satisfactory degree. Note that groundwater flow within the study area was assumed to be matrix flow, and local-scale conduit flow conditions were not considered, since the locations/dimensions of the subsurface cavities/voids, sinking streams, and springs were unknown due to lack of geophysical surveys.

Simulation Code
The MODFLOW computer code developed and released by the US Geological Survey was selected to develop the groundwater model in this study. MODFLOW was developed based on the concept of mass balance and Darcy's Law and is primarily utilized to simulate 3D constant-density groundwater flow through porous media (Harbaugh and McDonald, 1996). The governing equation is described in partial differential form: where, K xx , K yy , and K zz are values of hydraulic conductivity along the x, y, and z axes (L/T), h is the potentiometric head (L), W is a volumetric flux per unit volume representing sources and/or sinks of water (T −1 ), S s is the specific storage of the porous material (L −1 ), and t is time (T).

Spatial and Temporal Discretization
Spatially, the model domain was horizontally discretized into 248 rows and 218 columns with a uniform grid spacing of 30 m in both the x and y directions and was vertically divided into three layers, with Layer 1 representing the surficial aquifer, Layer 2 representing the upper confining unit, and Layer 3 representing the Floridan aquifer. The top elevations of Layers 1, 2, and 3 are shown in Figures 2A-C, and the bottom elevations of Layer 3 are shown in Figure 2D. Temporally, the model was steady-state under long-term annually averaged hydrologic and hydrogeologic conditions.

Parameters
For Layer 1, values of 30 m/day, 3 m/day, and 0.2 were assigned to horizontal hydraulic conductivity, vertical hydraulic conductivity, and porosity, respectively. For Layer 2, values of 0.01 m/day, 0.01 m/day, and 0.3 were assigned to horizontal hydraulic conductivity, vertical hydraulic conductivity, and porosity, respectively. For Layer 3, values of 600 m/day, 60 m/day, and 0.4 were assigned to horizontal hydraulic conductivity, vertical hydraulic conductivity, and porosity, respectively (Motz et al., 1995). Note that these values representing the characteristics of each layer were adjusted during the model calibration process.

Boundary Conditions
For Layer 1, general-head, no-flow, recharge, and evapotranspiration boundaries were assigned ( Figure 2E). The lateral boundaries that are parallel to the water table contours where groundwater flows into or out of the model domain were designated general-head boundaries. The lateral boundaries that are perpendicular to the water table contours where groundwater flux is zero were designated no-flow boundaries. The top of Layer 1 for representing infiltrated rainwater and groundwater evapotranspiration were designated recharge and evapotranspiration boundaries, respectively. Infiltrated rainwater was calculated by rainfall and infiltration/rainfall ratio (dependent on land use and land cover, which are shown in Figure 2G), and evapotranspiration rate was calculated by potential evapotranspiration rate.
Layer 2 was assigned a no-flow boundary. For Layer 3, general-head, no-flow, well, and drain boundaries were assigned (Figure 2F). The criteria for designating a boundary as general-head or no-flow were exactly the same as above and are not repeated here. A well boundary was assigned at abstraction wells, and a drain boundary was assigned at springs.

Initial Conditions
Initial heads were estimated based on the water table contours shown in Figure 1.

RESULTS AND INTERPRETATION Calibration
As mentioned above, the model was calibrated using the trial-and-error method by adjusting the values of the horizontal/vertical hydraulic conductivities within a reasonable range until the simulated groundwater levels and spring discharges matched the field measurements to a satisfactory degree ( Figure 3A). The field-measurements included: (1) the observation well (L-1020) that monitors potentiometric levels; (2) the observation station that monitors Mt. Plymouth Lake stages; (3) the observation stations that monitor the Droty Spring, Snail Spring, and Sulphur Spring discharges. Note that the time range of the observed groundwater levels and spring discharges is from 2010 to 2016, and the annual-average values are computed to serve as calibration targets since the temporal discretization of the groundwater flow model is steady-state.

DGLR and Relative Probability of Sinkhole Development
Based on the output from the calibrated model, the spatial variations of annual-average water tables, potentiometric levels, and DGLR are shown in Figures 3B-D, respectively. It can be observed that: (1) water tables are higher in the southwest and central portion and lower in the northeast and east portion (partly consistent with the topographic changes shown in Figure 2A); (2) potentiometric levels are higher in the southwest portion and lower in the northeast portion; (3) DGLR values are higher in the southwest and central portion and lower in the north, northeast, and east portion. Based on the simulated annualaverage DGLR, the relative probabilities of sinkhole development were identified as low (DGLR ≤ 50 mm/year), intermediate low (DGLR > 50 mm/year but ≤ 100 mm/year), intermediatehigh (DGLR > 100 mm/year but ≤ 150 mm/year), and high (DGLR > 150 mm/year), and spatial variation of relative probability of sinkhole development is shown in Figure 4A. Note that the locations of reported sinkholes (recorded by subsidence incident reports provided by the Florida Geological Survey) within the study area are also shown on Figure 4A (black dots). It can be observed that most of the reported sinkholes are located within areas where relative probability was classified as "high, " and the remaining ones are located within areas where relative probability was classified as "intermediate high, " indicating that the approach of using annual-average DGLR to identify the relative probability of sinkhole development in the study area is valid. The relative probabilities of sinkhole development at the under-construction expressway and its vicinity are shown in Figure 4B. It can be observed that more than half of the under-construction expressway (see the solid white line) is located within areas where the relative probability of sinkhole development was classified as "high" or "intermediate high, " indicating that engineering action (e.g., bridging over the "dangerous" areas) should be taken in advance to minimize the potential hazards of sinkhole development within the designed range of the pavement structures.

DISCUSSION
The simulated DGLR output from the developed/calibrated groundwater model can be used to generate a sinkhole susceptibility zonation map of the study area, and Figure 4B is a good example of such, serving as a useful indicator of the relative probability of sinkhole development. The development of regional-and local-scale groundwater models can contribute to the implementation of sinkhole warning systems to (1) help engineers and technicians to minimize the ground stability problems caused by sinkholes, (2) provide an information source for insurance companies and the public, and (3) implement a broad research platform for researchers.
There are some limitations to this study, along with a certain degree of uncertainty associated with simulation results. Firstly, this study ignored heterogeneity of the aquifer systems, which simplifies model implementation but sacrifices local-scale simulation accuracy. The spatial variations in hydrogeologic parameters are unknown due to lack of borehole data and pumping test data, and the aquifer systems are considered to be homogenous, with uniform values initially assigned to all model grid cells expected to be adjusted during calibration. Secondly, this study ignored local-scale conduit flow within the aquifer systems. The locations and dimensions of subsurface voids/conduits are unknown due to lack of geophysical surveys, and the real local-scale groundwater flow rate and directions may deviate from the simulated results. Thirdly, the spatial discretization is relatively "coarse" (30 m × 30 m). Finer spatial discretization would indeed improve computing accuracy and make model results more reliable and precise, especially as regards the hydraulic gradient at those areas where topographic change are huge and groundwater pumping rates are high. However, numerical instability is one of the problems that numerical simulations generally encounter. In order to lower the risk of numerical instability and maintain a reasonable computation runtime, one has to sacrifice computing accuracy to some extent. Fourthly, local-scale simulation results have uncertain range, since there is insufficient monitoring data on groundwater levels.
Note that evaluation of the relative probability of sinkhole development in this study is based on the finding that sinkholes are likely to occur in those areas that have higher DGLR. This finding has been demonstrated to be true in central Florida karst terrains in that it is developed based on local karst conditions, but it has not yet been verified in other places. In fact, the processes of sinkhole development are complicated and can be affected by multiple hydrological, geological, and geochemical factors such as heavy rainfall, rapid decline of groundwater level, thickness and composition of overburden sediments, surface loading, and subsurface transport of dissolved carbon dioxide. It is hard to tell whether DGLR would be the dominant factor affecting sinkhole development in other karst areas. To eliminate the geographical bias of this approach, future work would be extended to involve more impact factors to implement complex sinkhole development warning systems that can be applied for other vulnerable karst areas in west-central and north-central Florida.

SUMMARY
In this study, a groundwater model using the MODFLOW computer code was developed and calibrated against fieldmeasured groundwater levels and spring discharges to quantify annual-average DGLR for the purpose of evaluating the relative probability of sinkhole development at an under-construction expressway (built to connect two important state roads, SR46 and SR 429) and its vicinity in the central Florida sinkholeprone area. It was indicated that most of the expressway construction site has a relatively high probability of sinkhole development and that engineering action (e.g., bridging over the "dangerous" areas) should be taken in advance to minimize potential sinkhole hazards.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.