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 (Brinkmann et al., 2008; 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 over-exploitation, 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 annual-average DGLR. Xiao et al. (2016, 2018) conducted statistical analyses of the relationships between the spatial distribution of reported sinkholes and the spatial variation of annual-average 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.

FIGURE 1

Figure 1. Location of the study area and cross-sections A–A’ and B–B’ showing the hydrostratigraphic units (the bottom of the Floridan aquifer is much deeper than −50 m).

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²/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 field-measurements 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:

$$\frac{\partial }{\partial x}\left(K_{xx}\frac{\partial h}{\partial x}\right)+\frac{\partial }{\partial y}\left(K_{yy}\frac{\partial h}{\partial y}\right)+\frac{\partial }{\partial z}\left(K_{zz}\frac{\partial h}{\partial z}\right)+W=S_s\frac{\partial h}{\partial t}$$

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.

FIGURE 2

Figure 2. (A) Top elevations of Layer 1 (modified after US Geological Survey National Elevation Dataset); (B) top elevations of Layer 2 (modified after Williams and Kuniansky, 2016); (C) top elevations of Layer 3 (modified after Williams and Kuniansky, 2016); (D) bottom elevations of Layer 3 (modified after Williams and Kuniansky, 2016); (E) boundary conditions assigned to Layer 1; (F) boundary conditions assigned to Layer 3; (G) map of land use and land cover (modified from the GIS database maintained by St. Johns River Water Management District).

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.

FIGURE 3

Figure 3. (A) Locations of field-measurements (modified from the GIS database maintained by St. Johns River Water Management District) and scatter diagram showing the goodness of fit between the observed and simulated heads and spring discharges; (B) spatial variation of simulated water tables; (C) spatial variation of simulated potentiometric levels; (D) spatial variation of simulated DGLR.

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 annual-average 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), intermediate-high (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.

FIGURE 4

Figure 4. (A) Spatial variation of the relative probability of sinkhole development and locations of reported sinkholes; (B) relative probability of sinkhole development at the under-construction expressway (see the solid white line) and its vicinity.

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 field-measured 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 sinkhole-prone 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.

Author Contributions

Both authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

Funding

This research was partly supported by the Youth Program of the Natural Science Foundation of Tianjin (19JCQNJC08400), the Start-up Grant from Tianjin University of Science & Technology, and the Young Talent Fund of College of Marine and Environmental Sciences of TUST.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank Dr. Boo Hyun Nam and Dr. Dingbao Wang from the University of Central Florida (United States) for sharing their experiences and advice. We also thank the US Geological Survey, Florida Geological Survey, and St. Johns River Water Management District for providing data. The opinions, findings, and conclusions expressed in this manuscript were those of the authors and are not necessarily those of these agencies.

References

Anderson, M. P., Woessner, W. W., and Hunt, R. J. (2015). Applied Groundwater Modeling: Simulation of Flow and Advective Transport, 2nd Edn. Amsterdam: Elsevier.

Google Scholar

Beck, B. F. (1986). A generalized genetic framework for the development of sinkholes and Karst in Florida. USA. Environ. Geol. Water Sci. 8, 5–18. doi: 10.1007/bf02525554

CrossRef Full Text | Google Scholar

Brinkmann, R. (2013). Florida Sinkholes: Science and Policy. Gainesville, FL: University Press of Florida.

Google Scholar

Brinkmann, R., and Parise, M. (2009). The timing of sinkhole formation in Tampa and Orlando. Florida. Florida Geographer. 41, 22–38.

Google Scholar

Brinkmann, R., Parise, M., and Dye, D. (2008). Sinkhole distribution in a rapidly developing urban environment: hillsborough county, Tampa Bay area, Florida. Eng. Geol. 99, 169–184. doi: 10.1016/j.enggeo.2007.11.020

CrossRef Full Text | Google Scholar

Chitsazan, M., and Movahedian, A. (2015). Evaluation of artificial recharge on groundwater using MODFLOW model (case study: Gotvand Plain-Iran). J. Geosci. Environ. Prot. 3, 122–132. doi: 10.4236/gep.2015.35014

CrossRef Full Text | Google Scholar

Fleury, E. S., Carson, S., and Brinkmann, R. (2008). Testing reporting bias in the Florida Sinkhole Database: an analysis of sinkhole occurrences in the Tampa metropolitan statistical area. Southeast. Geogr. 48, 38–52. doi: 10.1353/sgo.0.0012

CrossRef Full Text | Google Scholar

Ford, D. C., and Williams, P. (2007). Karst Hydrogeology and Geomorphology. Chichester: Wiley.

Google Scholar

Gray, K. M. (2014). Central Florida sinkhole evaluation. Technical Publication, Florida Department of Transportation. District 5 Materials & Research.

Google Scholar

Gutiérrez, F., Parise, M., De Waele, J., and Jourde, H. (2014). A review on natural and human-induced geohazards and impacts in karst. Ear. Sci. Review 138, 61–88. doi: 10.1016/j.earscirev.2014.08.002

CrossRef Full Text | Google Scholar

Harbaugh, A. W., and McDonald, M. G. (1996). User’s Documentation for MODFLOW-96 – An Update to the U.S. Geological Survey Modular Finite-Difference Groundwater Flow Model, Open-File Report 96-485 (Reston, VA: USGS).

Google Scholar

Jammal, S. E. (1982). The Winter Park Sinkhole and Central Florida Sinkhole Type Subsidence. Winter Park, FL: Jammal and Associates, Inc.

Google Scholar

Kuniansky, E. L., Weary, D. J., and Kaufmann, J. E. (2015). The current status of mapping karst areas and availability of public sinkhole-risk resources in karst terrains of the United States. Hydrogeol. J. 24, 613–624. doi: 10.1007/s10040-015-1333-3

CrossRef Full Text | Google Scholar

Lindsey, B. D., Katz, B. G., Berndt, M. P., Ardis, A. F., and Skach, K. A. (2010). Relations between sinkhole density and anthropogenic contaminants in selected carbonate aquifers in the eastern Unites States. Environ. Earth Sci. 60, 1073–1090. doi: 10.1007/s12665-009-0252-9

CrossRef Full Text | Google Scholar

Mali, S. S., and Singh, D. K. (2016). Groundwater modeling for assessing the recharge potential and water tablebehaviour under varying levels of pumping and recharge. Indian J. Soil Conserv. 44, 93–102.

Google Scholar

Maroney, P. F., Cole, C. R., Corbett, R. B., Dumm, R. E., Eastman, K. L., Gatzlaff, K. M., et al. (2005). Final Report: Insurance Study of Sinkholes, Florida: Florida State Univ.

Google Scholar

Miller, J. A. (1986). “Hydrogeologic framework of the Floridan aquifer system in Florida and in parts of Georgia, Alabama, and South Carolina,” Paper Presented U.S. Geological Survey Professional Paper 1403-B. (Reston, VA: USGS).

Google Scholar

Motz, L. H., Beddow, I. I. W. D., Caprara, M. R., Gay, J. D., and Sheaffer, S. M. (1995). North-Central Florida Regional Ground-Water Investigation and Flow Model (Final Report), (Gainesville, FL: University of Florida).

Google Scholar

Sahoo, S., and Jha, M. K. (2017). Numerical groundwater-flow modeling to evaluate potential effects of pumping and recharge: implications for sustainable groundwater management in the Mahanadi delta region. India. Hydrogeol. J. 25, 2489–2511. doi: 10.1007/s10040-017-1610-4

CrossRef Full Text | Google Scholar

Sanford, W. (2002). Recharge and groundwater models: an overview. Hydrogeol. J. 10, 110–120. doi: 10.1007/s10040-001-0173-5

CrossRef Full Text | Google Scholar

Scanlon, B. R., Healy, R. W., and Cook, P. G. (2002). Choosing appropirate techniques for quanfitying groundwater recharge. Hydrogeol. J. 10, 18–39. doi: 10.1007/s10040-001-0176-2

CrossRef Full Text | Google Scholar

Tibbals, C. H. (1990). “Hydrology of the Floridan aquifer system in east-central Florida – regional aquifer system analysis – Floridan aquifer system,” Paper Presented U.S Geological Survey Professional Paper 1403-E. (Gainesville, FL: University of Florida)Google Scholar

Williams, L. J., and Kuniansky, E. L. (2016). “Revised hydrogeologic framework of the Floridan aquifer system in Florida and parts of Georgia, Alabama, and South Carolina” Paper Presented U.S. Geological Survey Professional Paper 1807. (Gainesville, FL: University of Florida).

Google Scholar

Wilson, W. L., and Beck, B. F. (1992). Hydrogeologic factors affecting new sinkhole development in the Orlando area. Florida. Groundwater 30, 918–930. doi: 10.1111/j.1745-6584.1992.tb01575.x

CrossRef Full Text | Google Scholar

Xiao, H., Kim, Y. J., Nam, B. H., and Wang, D. (2016). Investigation of the impacts of local-scale hydrogeologic conditions on sinkhole occurrence in East-Central Florida, USA. Environ. Earth Sci. 75, 1274–1289.

Google Scholar

Xiao, H., Li, H., and Tang, Y. (2018). Assessing the effects of rainfall, groundwater downward leakage, and groundwater head differences on the development of cover-collapse and cover-suffosion sinkholes in central Florida (USA). Sci. Total. Environ. 644, 274–286. doi: 10.1016/j.scitotenv.2018.06.273

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, Y., and Li, W. (2011). A review of regional groundwater flow modeling. Geosci. Front. 2, 205–214. doi: 10.1016/j.gsf.2011.03.003

CrossRef Full Text | Google Scholar