Assessment of groundwater quality and remediation in karst aquifers: A review

Author’s Accepted Manuscript Assessment of Groundwater Quality Remediation in Karst Aquifers: A Review

and

Koosha Kalhor, Reza Ghasemizadeh, Ljiljana Rajic, Akram Alshawabkeh www.elsevier.com/locate/gsd

PII: DOI: Reference:

S2352-801X(18)30163-2 https://doi.org/10.1016/j.gsd.2018.10.004 GSD163

To appear in: Groundwater for Sustainable Development Received date: 29 July 2018 Accepted date: 10 October 2018 Cite this article as: Koosha Kalhor, Reza Ghasemizadeh, Ljiljana Rajic and Akram Alshawabkeh, Assessment of Groundwater Quality and Remediation in Karst Aquifers: A Review, Groundwater for Sustainable Development, https://doi.org/10.1016/j.gsd.2018.10.004 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Assessment of Groundwater Quality and Remediation in Karst Aquifers: A Review Koosha Kalhor*, Reza Ghasemizadeh, Ljiljana Rajic, Akram Alshawabkeh Department of Civil and Environmental Engineering, Northeastern University, 360 Huntington Ave, Boston, MA 02115, USA *

Correspondence to: Koosha Kalhor , Current Address of First Author: 35 W35th St, 8th Floor, New York, NY 10001,[email protected]

Abstract Karst aquifers, capable of store and transmit large amount of water, are the main source of drinking water in many regions worldwide. Their excessive permeability leads to an enhanced vulnerability to store and convey contamination accordingly. From sustainability perspective, environmental, economic and social impacts are consequences of water pollution. In this study, an overview of hydrogeological processes and concepts regarding groundwater flow, and contaminant transport, are presented followed by a short discussion on surface water and groundwater interaction in karst aquifers. Due to the complexity of karstic systems, different approaches have been developed by researchers for investigating and predicting hydrogeological processes and groundwater behavior in karst which are reviewed herein. In addition, groundwater contamination and most common and effective remediation techniques in karstic terrains are discussed. In the last section, modeling techniques and remote sensing methods, as beneficial and powerful tools for assessing groundwater flow and contaminant transport in karst terrains, are reviewed and evaluated. Moreover, in each section, associated research works conducted for Puerto Rico are discussed and some recommendations are presented to complement ongoing hydrogeological investigations on the island. Graphical Abstract:

Keywords: Groundwater; Karst; Remediation; Hydrogeology; Modeling; Puerto Rico

1 Introduction Sustainable water resources management is a crucial concern in most countries across the globe. Only 3 % of total water on the Earth is considered as fresh water resources and approximately 30 % of that is accessible as groundwater, which is vital for human health, ecosystem, energy industry and other water-dependent topics (Shiklomanov, 1993). Karst aquifers are responsible for providing potable water for 40 % and 25 % of the US and world’s population, respectively (Ghasemizadeh et al., 2012). The increasing demand by residential, industrial and agricultural uses have caused groundwater depletion and decreased quality in many regions. From sustainable development perspective, environmental, economic and social impacts are consequences of water pollution in any region of the world. Hence, careful attention should be paid to preserve water resources. With particular reference to karst aquifers, this study provides a comprehensive review of hydrological concepts and novel investigation and modeling techniques followed by a discussion of groundwater contamination and remediation

techniques. This paper aims to present and review the work of other researchers in the recent years (especially after 2010) and to discuss the latest improvements that have been made regarding groundwater quality and quantity assessment. In each section, the associated research work that has been done for Puerto Rico is presented to better understand what research works have been conducted and what are the research gaps on karst aquifers of this island. Puerto Rico, as the case study location, is considered a territory of the United States (US). The island is located in northeastern side of Caribbean Sea and have an estimated population of 3.6 million (Castro-Prieto et al., 2017). Several surface water and groundwater resources across the island provide residents with fresh water and are used for agricultural, industrial and energy-based purposes. Figure 1 exhibits the geographical location of Puerto Rico, generalized geological setting of the island, surficial geological formation of north coast limestone aquifer and crosssectional view of geological layers along the north coast.

Figure 1. Geographical location of Puerto Rico (a), generalized geological setting of the island (b), surficial geological setting of north coast limestone aquifer (c) and cross section view of geological layers in Florida and Barceloneta municipalities (Vertical scale is exaggerated, d) – modified from (Gómez-Gómez et al., 2014)

Due to the presence of karstic aquifers with high level of heterogeneity and anisotropy in northern coast of the island, water from precipitation can rapidly percolate into underground. This rapid movement, makes karst aquifer highly vulnerable to emerging contamination from several sources of contamination induced by agricultural and industrial activities in addition to proximity to urban areas (Cherry, 2001). Moreover, highly heterogeneous karstic aquifers with a lot of conduits cause high rate of water level fluctuation even in small temporal scaling. Yu et al. assessed the patterns of temporal scaling of groundwater level fluctuation for the karstic aquifer of Puerto Rico and pointed out that according to multifractal or singularity spectrums, there are smoother fluctuations in the alluvial aquifer and rougher ones in the northern karst aquifer (Yu et al., 2016).

2 Karst Aquifers 2.1 Characterization Comprising of chemically soluble rocks with large passages or network of conduits and caves inside, karst aquifers are very permeable and capable of storing and transmitting large amount of water. Therefore, karst regions have great potential of providing habitants with fresh water (Quinn et al., 2006). In fact, 20-25 % of the world’s population depends on water supplies from karst aquifer directly or indirectly. Approximately 10 % of the world’s land surface areas have karst aquifer beneath them. This percentage is higher in some areas such as in Europe where it is roughly 35 % (Ford and Williams, 2007). Until recent years, the boundaries of karst aquifers around the world were not recognized accurately. Hence, by taking advantage of GIS tools, recent exploration of karst aquifers and Global Lithological Map that was developed before, Chen et al. completed the first World Karst Aquifer Map (WOKAM). Their map distinguishes continuous carbonate rocks and discontinuous carbonate rocks and include major karst springs, wells and caves (Chen et al., 2017). Many other resources show karst map of the US or the world. Distribution of karst aquifers within the United States and its territories based on USGS data is presented herein (Fig. 2). Detailed geological investigations in Puerto Rico show that limestone, dolomite, gypsum and anhydrite are the most common materials forming Puerto Rican karst aquifers (Fig. 1).

Legend Carbonate Rocks Evaporite Rocks Sedimentary Rocks Quartz Sandstone Volcanic Rocks Evaporite Basins

Figure 2. Distribution of karst aquifers in the United States and its territories – Compiled from open files associated with the USGS report of (Weary and Doctor, 2014) Karst aquifers are individually different with unique work-frame and characteristics and they should be studied case by case (Stevanović, 2015). Two important characteristics of karst aquifers are heterogeneity and anisotropy which make it hard for hydrogeologists and researchers to develop models using simplifying assumptions. Karsts have the most complex system amongst terrains, causing lot of uncertainties and errors in developed models for studying and predicting their behavior (Bakalowicz, 2005). Also, the recharge and discharge rate of karst springs can vary a lot due to several reasons such as fluctuations in water table level caused by hydrological events or seasonal variations (Gárfias-Soliz et al., 2009). Table 1 elaborates hydrogeological characteristics of three main aquifer types, porous media, fractured rock and karst system based on ASTM D 5717–95 Standard (ASTM, 1995).

Table 1. Hydrogeological characteristics of three main aquifer types, porous media, fractured rock and karst system based on ASTM D 5717–95 Standard (ASTM, 1995; Rosenberry and LaBaugh, 2008) Aquifer characteristics Effective porosity

Porous (Granular) Mostly primary, through intergranular pores

Isotropy Homogeneity Flow

More isotropic More homogeneous Slow, laminar

Flow predictions Storage

Darcy's law usually applies Within saturated zone

Aquifer type Fractured Rock Mostly secondary, through joints, fractures, and bedding plane partings Probably anisotropic Less homogeneous Possibly rapid and possibly turbulent Darcy's law may not apply Within saturated zone

Karst Mostly tertiary (secondary porosity modified by dissolution); through pores, bedding planes, fractures, conduits, and caves Highly anisotropic Heterogeneous Likely rapid and turbulent Darcy's law rarely applies Within both saturated zone and epikarst

Recharge

Dispersed

Primarily dispersed, with some point recharge

Temporal head variation Temporal water chemistry variation

Low variation Low variation

Moderate variation Low to moderate variation

Ranges from almost completely dispersed- to almost completely pointrecharge Moderate to high variation Moderate to high variation

A karst aquifer system comprises several elements such as caves, conduits, sinkholes and springs. Limestone karst aquifers are common in many areas around the world including Puerto Rico (Cherry, 2001; Rafael et al., 2016), Florida (Xu et al., 2016) Mexico (Bauer-Gottwein et al., 2011) and China (Luo et al., 2016). Karsts usually are evolved from fractured or fractured-porous rock networks after several years through carbonate dissolution creating large passages and caves. Several modeling methods have been employed to simulate the evolution of karst aquifers from fractured or porous-fractured rock systems (Kaufmann, 2003, 2016; Kiraly, 2003). 2.2 Hydraulic Conductivity In groundwater hydrology, hydraulic conductivity quantifies the ability of soil in transferring water. Based upon various aquifer materials, hydraulic conductivity can range from 10 cm/s for gravel to 10-10 cm s-1 for shale. Hydraulic conductivity of karst limestone is the greatest comparing to many other aquifer types (Fig. 3).

Karst limestone Permeable basalt

Fractured metaphoric and igneous rocks Limestone and dolomite Sandstone Unfractured metaphoric and igneous rocks -11

-10

-9

-8

-7

-6

-5

-4

-3

-2

-1

0

1

Log K (cm/s) Figure 3. Hydraulic conductivity (K) range for different types of rock – modified from (Freeze and Cherry, 1979) Laboratory, field and numerical methods are 3 main means of measuring hydraulic conductivity. Numerical and finite element-based methods are used for determining vertical and horizontal hydraulic conductivities (Kalbus et al., 2006; Smith et al., 2016). Usually, in karst aquifers where subsurface heterogeneity exists, determining hydraulic parameters such as K requires a complicated analysis because this parameter is spatially and temporally variable throughout the aquifer. Hydraulic tomography is a novel method that can be used for imaging the heterogeneity in karstic terrains (Illman et al., 2007). Moreover, in karst aquifers, the average range of K can vary depending on several factors such as geology, slope, level of heterogeneity and karstification. Angulu et al. studied hydraulic conductivity in karst areas by applying water injection tests and electrical resistivity logging (Angulo et al., 2011). Similarly, different researchers reported experimental values for hydraulic conductivity based on their research approach and case study area (Chen et al., 2011; Fu et al., 2015; Sudicky et al., 2010). As it is shown in Fig. 3, K of limestone karst aquifer which is dominant in northern coast of Puerto Rico, can be assumed in the range of 10 -4 to 5 cm s-1 which demonstrates high level of

permeability in karst aquifers. The estimated values of K in different locations of north coast karst aquifer of Puerto Rico can be found in the water resources investigation reports (Rodriguez-Martinez, 1995) or similar sources for fine-scale studies or modeling purposes. 2.3 What methods can be used to study karst aquifers? Based on the aforementioned complex characteristics of karst, several techniques and tools associated with modified and reformed conventional methods (such as hydrologic and hydraulic methods, geophysical and geological methods, modeling techniques and tracer tests) have been employed by researchers for understanding the behavior of karst aquifers. (Goldscheider and Drew, 2007; Stevanović, 2015). (Giudici et al., 2012; Hu et al., 2009) studied karst aquifers by taking advantage of modeling methods. In addition, taking advantage of geological methods, which help understanding the aquifer geometry and hydraulic properties such as permeability in addition to orientation and characteristics of potential flow paths, can boost the accuracy of modeling results (Goldscheider and Drew, 2007). Geophysical techniques can also be employed in conjunction with geological methods to understand geologic structures and overburden thickness of the aquifer (Chalikakis et al., 2011; Ford and Williams, 2007; Goldscheider and Drew, 2007). 2.4 Surface Water and Groundwater Interaction (SWGWI) in Karst Interaction of surface water (SW) and groundwater (GW) plays a critical role in understanding hydrological behavior of a basin. This interconnectivity incorporates the topographical, geological and morphological characteristics of terrains. Generally, water recharge from inflow of GW into the riverbed, water discharge from river bed to aquifer and also losing and gaining water for both SW and GW in some river segments are three main processes that can occur in SWGWI (Winter et al., 1998). High permeability and low attenuation capacity of karst aquifers make mixture of SW and GW problematic for fresh water use in karst terrains. Rainfall infiltration and streamflow response to base-flow discharge are separated by a much shorter delay in karst compared to other aquifer types (Bayless et al., 2014). Due to the complex patterns of SW and GW flow in karst, several studies have denied coincide of SW and GW drainage divides. A river that disappears in a SW basin and reappears in another basin can be mentioned as an example which depicts uncertainties in understanding the source of water and its associated dissolved contaminants in karst (Winter et al., 1998). Usually, carbonate rocks are present at the land surface in karstic terrains; however, in some cases, the bedrock is covered by other deposits (e.g. surficial deposits) which is called “mantled” karst (e.g., Edwards Aquifer in south-central Texas). In mantled karst, shallow GW has interaction with lakes and wetlands in a way similar to that in sandy glacial and dune terrains. The small difference in interaction of these two can be tied with the buried carbonate rocks. Dissolution of buried carbonate rocks can lead to slumpage of an overlying confining bed, allowing water to flow readily through the confining bed. Consequently, changes in hydraulic heads within the aquifers underlying the confining bed can affect the lakes and wetlands (Winter et al., 1998). In addition to south-central Texas, mantled karst can be found in north-central Florida. This area has several sinkhole lakes and unconsolidated deposits, overlie the highly soluble limestone of the Upper Floridan aquifer. When unconsolidated surficial deposits slump into sinkholes which form in the underlying limestone, most land surface depressions containing lakes are made. Hence, in most of the times, sinkholes in the bedrock underlying lakes affect the hydrology of the lakes, although the lakes are not located directly within limestone (Winter et al., 1998).

The general term “Surface water” can be referred to wastewater on the ground, brackish water of the sea/ocean, and water in lakes and rivers. Attempts have been made to minimize the interaction between polluted SW and GW in karstic areas by placing an impermeable seal along canal bottoms or riverbeds, constructing small dams/weirs, building grouting curtains, changing the flow direction of surface waters, plunging ponors and creating reactive barriers by pumping freshwater into aquifer (Milanović et al., 2015). Modeling, hydrographic analysis and hydrogeological mapping are among common “desktop” techniques which are used to assess SWGWI in karst (Fleckenstein et al., 2010). Hydrographic analysis, which monitors stream flow time series and define base flow component, can be used easier compared to the other two methods. However, it is only applicable for gaining stream conditions. Hence, it may not give accurate result if used for a karst aquifer. Hydrogeological mapping method, which provides comprehensive understanding of GW systems around streams and their related hydrogeological systems, can describe hydrogeological components of a GW system at a course scale (Gleeson et al., 2014). Thus, it can be used for regional assessment of karst aquifers. In addition, modeling methods can simulate water flow and contaminant transport around streams and within aquifer using mathematical equations. Although they can predict the behavior of a hydrogeological system, their application in karst is limited. Usually, simplifying assumptions are made during modeling process of SWGWI in karst which result in uncertainties. However, by taking advantage of new techniques/technologies and/or by combing this method with other in-situ measurements, researchers have shown the application and promising future of modeling methods (Guay et al., 2013).

Also, field observations, seepage measurement, hydrochemical methods such as environmental isotopes, hydrochemistry, tracers and also hydrometric and geophysical analysis are beneficial “field” tools which can be used for describing SWGWI in karst (González-Pinzón et al., 2015; Martinez et al., 2015). In the recent years, use of methods associated with geophysics (e.g. resistivity, EM, radiometrics) or remote sensing (e.g. Landsat) has been reported for mapping landscape features that indicate or control connectivity (Meyerhoff et al., 2013). However, because these methods require specific equipment, technical expertise and logistical support, researchers have had tendency to use more practical and “easy to use” methods such as Hydrometrics and Water Budgets analysis (Brodie et al., 2007). Sometimes, because of unknown and complex catchment boundaries, analyzing water budgets can also be problematic (Goldscheider and Drew, 2007; Hartmann et al., 2014; Kovács et al., 2005). Among field tools, tracers have been widely used due to their capability of providing independent ways of validating or refuting conventional-traditional methods of analyzing data and describing SWGWI (Baskaran et al., 2009; Jankowski, 2007). By using tracers such as florescent dyes, aquifer parameters such as recharge and discharge and fluid transport properties can be determined (Martinez et al., 2015). However, tracer studies require careful planning (e.g. meeting environmental regulatory controls) and dosage control. In many cases, isotropic techniques and artificial tracers are used for determining residence time and water age and for understanding the water movement through conduits. The main advantages of these techniques are determining linear flow velocities and information on contaminant transport and delineating catchment areas. Although obtained information and data from tracers are often reliable and unequivocal, limited applicability in large areas with long transit times and also change of color and toxicity concerns are some of the disadvantages of using isotropic techniques and artificial tracers (Goldscheider et al., 2008; Jones and Banner, 2003; Morales et al., 2017). Reviewing the literature, it seems like the use of tracers has been more frequent and beneficial for understanding the interconnectivity between SW and GW in karst terrains compared to other field methods (Katz et al., 1997), for example evaluating the connection

between sinking streams and down-gradient springs. It should also be noted that temperature change analysis and water budget assessment can be coupled with other methods to achieve more accurate and validate results (Brodie et al., 2007). Due to complexity of karst’s behavior, different types of “integrated” techniques and/or less common approaches have been employed by researchers to evaluate their performance in describing the interconnectivity between SW and karst GW (Brian Katz et al., 1997; Chu et al., 2016; Loáiciga, 2001; Rugel et al., 2016). Sutton and Screaton described SWGWI of a karst aquifer basin in Florida by analyzing river discharge and a transient numerical GW flow modeling. Their modeling results highlight the prominence of spatiotemporal variations in head gradients that can affect streams and karst aquifers connections and aquifer martial dissolution (Sutton et al., 2014). Bailly-Comte et al. studied the hydrodynamic interactions between GW and surface water in a karst watershed in southern France. Their focus was on the effect of GW on the genesis and propagation of surface floods. They showed the role of initial water level in a karst aquifer in predicting the type of hydraulic connection between SW and GW during flood events (Bailly-Comte et al., 2009). Moreover, the analysis of SW and karst GW interaction conducted by (Bayless et al., 2014) showed that analytical methods such as hydrograph separation and hysteretic loops can be used for identifying bounding conditions within the watershed. Chu et al. studied SWGWI in a karst aquifer in China using GW levels data and hydrochemistry analyses, together with isotope data based on hydrogeological field investigations. Their cross-correlation analysis showed that precipitation results in increasing GW level with a 2-month lag time. Moreover, their spectral analysis depicted the possible inter-connectivity between GW level, GW exploitation and precipitation (Chu et al., 2016). Rugel et al. conducted longitudinal sampling on a creek in Georgia, USA. They realized there are stepwise GW inputs from the Upper Floridan Aquifer, (in addition to specific conductance increase) at 5 out of 50 sampled reaches. The results of sampling along the 50-km study river section revealed that 42 % of total GW inputs entered through discrete preferential flow pathways were located within these 5 sampling locations (Rugel et al., 2016). In north coast limestone aquifer of Puerto Rico, with presence of rivers and lagoons and occurring intense precipitation and seawater intrusion, SWGWI can have a major impact in the quality of fresh water within karst aquifers of the region. However, we failed in our attempts to find any detailed research works associated with SWGWI in karst aquifers of northern PR. Hence, collecting field data, doing research and developing models to fully understand the mechanism of SWGWI in northern part of the island is strongly recommended. Despite the fact that numerous methods have been developed for describing SW-GW interrelationship, there are still uncertainties and lack of sufficient knowledge for fully understanding 1- the time lag between GW pumping and its influence on SW, 2- relationship between GW pumping and river losses and 3- exact recharge and discharge points in streams (Jankowski 2007). Wu et al. assessed uncertainties in SWGWI modeling. Employing a probabilistic collection method, they evaluated the applicability of the frame-work through an integrated SW-GW model for a basin in China and asserted that in describing complex SWGWIs, modeling uncertainties depend on the output and have significant spatiotemporal variability. Hence, employing a systematic uncertainty analysis can be extremely helpful in understanding SWGWI (Wu et al., 2014) and also in groundwater model development process of karst aquifers (Engelhardt et al., 2014). 3 Groundwater Contamination and Remediation Techniques

Karst aquifers which are portrayed as high permeable soil/rock systems with caves and fractures inside and also recharged by sinkholes and rivers, have shown high vulnerability to contamination (Kačaroğlu, 1999). The formation of solution channels and sinkholes facilitates the intrusion of seawater and contaminated storm water and wastewater into the aquifer. This is the main reason of accelerating the spread of contaminants in karst more than other aquifer types (i.e. porous media and fractured rock). Basically, fate and transport processes are not only different for each contaminant (depending on the physical and chemical properties of the contaminant), but also are particularly different in karst compared to other aquifers. A few studies assessed historical GW contamination in north coast limestone karst aquifer of Puerto Rico and suggested beneficial remedial actions and water resources management approaches (Biaggi, 1995; Padilla et al., 2011). 3.1 Contaminant types and sources of contamination Agricultural, industrial, residential, commercial and municipal development are considered as the main sources of groundwater pollution in recent decades (Fetter, 2001; Wakida and Lerner, 2005). In particular, leakage of storage tanks, chemical spills, landfills, fertilizers and pesticides, sanitation systems, untreated waste discharge and sewage are some of the main sources of contamination due to anthropogenic activities (El Alfy and Faraj, 2017). Generally, regardless of the source of contamination, pathogens, organic compounds (Lapworth et al., 2012), metals (Yao et al., 2012), and other inorganic compounds (e.g., nitrates, chlorides) are often found in groundwater (Panagiotakis and Dermatas, 2017; Vidal Montes et al., 2016). The most commonly found contaminants in karst GW are presented below. 3.1.1 Nitrate Nitrate is one of the most common GW contaminants according to several studies (Almasri and Kaluarachchi, 2004; Babiker, 2004). Point sources of Nitrogen/Nitrate leachate such as old septic systems, landfills, wastewater holding ponds and leaks from cracks in sewer pipelines can cause GW Nitrate contamination (Almasri, 2007; Kendall and Aravena, 2000; Wakida and Lerner, 2005). Additionally, non-point sources of N leaching, such as fertilizer use in agricultural areas, play a very significant role in increasing GW Nitrate contamination. In a USGS report with focus on Manati and Vega Baja municipalities in Puerto Rico, Nitrate occurrence and contamination were assessed. It was identified that the major sources of Nitrate contamination in the karst aquifer of the region are use of fertilizers for cultivation of pineapples and also septic tank effluent in rural and un-sewered (no sewer system) areas (Conde-Costas and Gómez-Gómez, 1999). 3.1.2 PPCP Nowadays, treating water contaminated by new chemical compounds mainly originated from pharmaceutical and personal care products (PPCP) are a big concern (Metcalfe et al., 2011). Concentration of PPCPs such as Antibiotics, Anti-inflammatories, Lipid regulators, Psychiatric drugs, Stimulants, Insect Repellants and Sunscreen agents is higher than regulatory criterion in many areas. Usually, wastewater and contaminated surface water, landfills, septic systems and sewer leakages are considered as common sources of PPCP contamination especially in karstic areas (Dodgen et al., 2017; Sui et al., 2015). High concentration of PPCPs in karst aquifers of Switzerland and Jordan was reported by (Morasch, 2013) and (Zemann et al., 2015), respectively. Moreover, Dodgan et al. assessed PPCP concentration in karst aquifers of southwestern Illinois in the USA and found out that no single water quality or meteorological parameter was significantly associated with PPCP contamination across all site profiles. They mentioned septic tanks effluent as a probable source of PPCPs (Dodgen et al., 2017). In another study, it was realized that due to high transmissivity of the

shallow alluvium, the GW below alluvial deposits is subjected to pollution caused by PPCP residues from wastewater effluents through artificial aquifer recharge (Sui et al., 2015). No studies related to PPCP contamination of GW in PR were found. However, it is expected that due to large number of septic tanks and poor effectiveness of wastewater treatment plants to deal with PPCPs in addition to use of mosquito repellent creams (because of Zika virus), groundwater quality will be affected by PPCP contamination in the island (Hennessey et al., 2016). 3.1.3 Metals In carbonated water of the karst aquifer (at natural pH), heavy metals such as Cr, Ni, Pb and Cd often precipitate as hydroxides and carbonates. Heavy metal transport, as an episodic phenomenon, mostly occurs when metals are adsorbed onto small particles of soil. When hydraulic condition of GW flow facilitates suspension of small particles, heavy metals can migrate down the flow path (Vesper et al., 2003). Also, elevated concentration of metalloids such as Arsenic has raised concerns in many regions. Heavy metals and metalloids in groundwater can be present in areas where there are mining operations or mine waste dumps. Industrial activities and urban waste are known as possible sources of metal(loid)s leachate in karstic GW (Alloway, 2013). Additionally, urban surface runoff containing high concentration of metals goes through karst aquifers via sinkholes and conduit network (Bonneau et al., 2017). Heavy metal ions are also frequently leached naturally from rocks and soils within karst media and can be introduced with acidic deposition (Drew and Hötzl, 1999). In karst GW of Puerto Rico, high concentration of metals (mainly arsenic, lead, and total chromium) have also been detected mainly at superfund sites (Fig. 4) (USEPA, 2003). 3.1.4 LNAPL and DNAPL Light Non-Aqueous Phase Liquids (LNAPLs) have tendency to float on the surface of groundwater table. They are usually present behind obstructions in cave streams within karst aquifers. Unlike LNAPLs, Dense Non-Aqueous Phase Liquids (DNAPLs), sink to the bottom of the aquifer and have tendency to be stored in lower points of the bottom of conduits. Moreover, it was reported that they tend to concentrate in open conduits (Ghasemizadeh et al., 2012). LNAPL and DNAPL Transport in karst aquifers is contingent upon runoff and GW flow hydraulics. LNAPLs can be forced into the karst system as plug flow and DNAPLs can be mobilized by transporting sediments at the bottom layer. In contrast to dissolved contaminants, the transport of NAPLs in aquifers is strongly influenced by buoyancy effects and surface tension properties. This makes the transport behavior of NAPLs less directly related to groundwater flow and more difficult to predict (Vesper et al., 2003). By conducting transport experiments using non-aqueous phase TCE in a laboratory-scale karst limestone rock, Carmona De Jesús and Padilla showed that although DNAPLs prefer to move downward, they can also move along non-vertical conduits when subjected to GW flow fields. Furthermore, it was demonstrated that once DNAPLs penetrate preferential flow pathways, subsequent injections often follow the same route. When they reach limited flow areas, DNAPLs will be accumulated and transported through diffusive processes after dissolution (Carmona De Jesús and Padilla., 2015). 3.1.5 Volatile Organic Compounds (VOC) Volatile and Semi-Volatile Organic Compounds (VOCs and SVOCs, respectively) flow along with the subsurface waters in karst; but their volatility, results in migrating the vapor phases of these contaminants through the vadose zone. (Field, 2017). A few studies assessed the contamination of Chlorinated Volatile Organic Compounds (CVOC) in PR (Padilla et al., 2011; Rivera et al., 2017). Historical studies since 1980 show that mainly, contaminants with chlorinated solvents

including TCE, dichloroethene, chloroform, carbon tetrachloride, tetrachloroethene, tetrachloroethane, dichloroethane and methylene chloride, have high concentrations causing public health concerns (Padilla et al., 2011). Several wells and sites (Fig. 4) were considered as the National Priority List (NPL) Superfund Sites and remediation actions for treating water in these sites begun years ago (Yu et al., 2015).

Regarding studying the groundwater pollution and understanding the potential exposure pathways of contaminants, some methods such as using tracers and GIS has been employed by researchers (Steele-Valentín and Padilla, 2009). The abundancy of superfund sites and high concentration of contaminants in GW recourses of Puerto Rico are suspected to be responsible for increasing rate of pre-term birth (highest amongst US states and territories) in the island (Mathews and MacDorman, 2011). However, since 2006, this rate has been declined from 20 % to 11.4 % due to remediation actions and enhanced awareness of habitant regarding water-borne diseases (March of Dimes Website, 2016; Rutigliano, 2016).

Figure 4. Case study location (upper picture), RCRA, NPL Superfund sites, aquifers, and the sampling wells (lower picture) based on the study of Yu et al. 2015 As it appears in Fig. 4, there is abundance of superfund sites in NPR mainly due to industrial activities, improper management of landfills, accidental spills, unidentified waste disposals, or residential septic systems. Most of these sites are located in upper aquifer of North Coast Limestone aquifer system. On the border of Arecibo and Barceloneta, there are 3 superfund sites which indicated high level of contamination in groundwater of that area. In karst aquifers of PR, analysis of contaminant transport requires more sophisticated approaches coupled with continuous field assessments and data collection (Hu et al., 2009). Yu et al. assessed the concentration of CVOCs such as Trichloroethylene (TCE), Perchloroethene (PCE), Carbon Tetrachloride (CT), Chloroform or Trichloromethane (TCM), and Methylene Chloride (DCM), in northern Puerto Rico based on historical data. They stated that the hydrogeological

conditions of the karst aquifer (intrinsic properties and the biological environment) in addition to the source origin were greatly associated with the spatiotemporal distribution patterns of the CVOCs. Since water resources pollution in northern Puerto Rico has caused a negative social, economic and environmental impact, long-term, consistent monitoring of water quality in addition to implementation of remedial actions were suggested (Yu et al., 2015).

3.1.6 Pathogens Pathogens such as viruses and bacteria are easily transported in karstic systems. This is because of high drainage potential of karst and lack of soil layers, which can act as filters. Hence, pathogens retain their activity in karst conduits for a long time (Vesper et al., 2003). Fecal coliforms, fecal streptococci, and often a particular emphasis on Escherichia coli have been studied more compared to viruses in karst aquifers. These organisms often indicate GW contamination caused by sewage or animal waste. Studies regarding bacterial contamination in West Virginia (Mathes, 2000), Ontario (Conboy and Goss, 2000), and western Ireland (Thorn and Coxon, 1992), show the presence of fecal bacteria in carbonate wells (White, 2016). Muldoon et al. assessed sources of fecal contamination in the Silurian dolomite aquifer of Wisconsin using enteric pathogens (Muldoon et al., 2016). Moreover, it was found out that high percentage of bacteria travels within karst media by adsorbing on small particulates (mostly clay) (Mahler et al., 2000; White, 2016). Resources for assessment of pathogens in karst aquifers of Puerto Rico are scarce. By assessing the North Karst Belt Zone of Puerto Rico, Acosta-Colón suggested that mesofauna diversity can be a possible indicator of pathogenic and opportunistic species. It was also recommended that in addition to characterizing mesofauna and species, more knowledge about the diet of mesofauna organisms is needed to clearly identify species that can be used as indicators of cave microbial pathogens (Acosta-Colón, 2016).

3.1.7 Chloride Coastal aquifers, such as north coast limestone aquifer of Puerto Rico, are susceptible to seawater intrusion which can increase the salinity (Chloride concentration) of groundwater (Arfib et al., 2007; Mongelli et al., 2013). Having a hydraulic connection to the sea, karstic-coastal aquifers can be characterized by having groundwater flow in conduits, sub-marine freshwater springs and intrusion of seawater through the aquifer via conduit networks (Fleury et al., 2007). Due to susceptibility of karstic-coastal aquifers, developing a sustainable plan by using integrated models for managing and monitoring water resources within these aquifers is essential (Sreekanth and Datta, 2015). Same as pathogens, not much research has been conducted regarding seawater intrusion in coastal karst aquifer of Puerto Rico. Although Zack et al. stated that the freshwater-saltwater interface of the North Coast limestone aquifer of PR is probably seaward and freshwater discharges to the seabed, by extracting water from public supply wells which results in decreasing GW level, seawater intrusion can be a potential problem in NPR (Zack et al., 1987). Not much chloride concentration data is available within USGS database and other resources. Hence, conducting a field sampling for a few years is recommended to assess the potential risk of chloride-based GW contamination. 3.2 Remediation strategies in karst aquifers

Based on increased understanding of contaminant transport behaviors and efficacies of various remediation technologies, USEPA emphasizes that combination of source and dissolved phase (plume) remediation is preferable strategy for remediation of many sites (USEPA, 2013). However, in addition to the challenges in costs and the extend of the treatment efficiencies arising at different sites, remediation in karst is even more challenging than in other

Percentage of in-situ source treatment decision documents

hydrogological environments due to the presence of highly permeable conduits of unknown extent (Parise et al., 2015).

60

(a) 50 40 30 20 10 0 2005

2006

2007

Percentage of GW Remediation Decision Documents

100

2008 Year

2009

2010

2011

(b)

90 80 70 60 50 40 30 20 10 0 1985

1988

1991

1994

1997

2000

2003

2006

2009

2012

Year Containment (VEB) In-situ Treatment Pump & Treat

MNA Other GW Remedies

Figure 5. Preferred remediation techniques for a) source zone treatment and b) dissolved phase (plume) treatment implemented at superfund sites (MNA and VEB stand for Monitored Natural Attenuation and Vertical Engineered Barrier, respectively. Other GW Remedies include Institutional control) – modified from (USEPA, 2013) For example, Sinreich stated that at the field scale, selecting the specific remediation technique for each contaminant requires comparative tracing experiments using contaminant surrogates (Sinreich, 2014). In spite of the dominant preferential flow pathways, solute and colloid tracers interact with aquifer material. Hence, if reactive and/or nonpersistent contaminants are present in conduits, hydro-chemical properties of contaminants are more important than aquifer intrinsic vulnerability in the transport process and arrival of contaminants at karst springs. Most remediation

processes in karst have lower efficiency compared to other aquifer types because of their time-dependent characteristics. Sinreich also suggested that the assumption of high vulnerability of karst may not be true for all conditions and should be identified based on fate and transport of specific contaminants. (Sinreich, 2014). Thus, no preference should be considered for karst GW treatment methods because they can have different efficiency and cost based on contaminant properties and site condition. Another example related to pump and treat application (Field, 2017), confirmed varying levels of success due to extreme heterogeneity and anisotropy of karst, limited diameter/length of the conduits and the behavior of flow regime within the solution conduits. It was realized that pump and treat was barely used by agencies for remediating contaminated karst GW (USEPA, 2013). As concluded by Parise et al., main difficulties with achieving feasible plume containment in karst with pump and treat systems are finding all preferential flow paths of the dissolved contaminant to the points of compliance and hydraulically stopping this transport by pumping both contaminated and clean water (Parise et al., 2015). Due to all these challenges, it was pointed out that contaminated karst aquifers tend to fall under the special category of Technical Impracticability (TI) waivers. This means regardless of how much effort and financial resources are expended, actual remediation cannot be accomplished because of technical implausibility especially for aquifers contaminated by DNAPLs (Field, 2017; USEPA, 2012). It was found that most of the case studies associated with TI waivers were those in karstic hydrogeological setting (Deeb et al., 2011).

However, in situ thermal treatment, in situ chemical oxidation, and in situ bioremediation are increasingly being applied at fractured rock and karst site for source and plume remediation while point-of-use and point-of-entry has been evaluated and proposed strategy for karst groundwater treatment (Randrianarivelo et al., 2017). 3.2.1 Source and Plume Remediation In-situ thermal treatment, In Situ Chemical Oxidation (ISCO), and in situ bioremediation are the main in situ remediation methods which have commonly been used in karst aquifers with certain limitations (Beyke et al., 2014; Parise et al., 2015). Studies have shown that identifying the location of conduits and major fractures is necessary for an efficient remedial action. It was reported that Electrical Resistance Heating technique was successfully implemented for removing DNAPLs and TCE in a karst aquifer in Alabama but preferential pathways within karst aquifers leading to high seepage velocity can cause heat loss (Hodges et al., 2014). The preferential flow pathways within karst system which can disperse injected materials, usually are problematic for other treatment methods such as ISCO and bioremediation as well.

Regarding in situ bioremediation, Hashim et al. pointed out that biological or biochemical methods which use microbes and nutrients for bio-precipitation, enzymatic oxidations, bio-surfactants and sulfate reductions for removing heavy metals, seems to be applicable and efficient in karst because injecting nutrients and electron donors is relatively inexpensive and non-toxic (Hashim et al., 2011). Byl and Williams stated that in areas isolated from the major groundwater flow paths, low migration of contaminant may possible to allow biodegradation of chlorinated ethenes and fuel if favorable microorganisms, food sources, and geochemical conditions are present (Byl et al., 2001; Byl and Williams, 2000). The potential for biodegradation of chlorinated organic solvents was in a karst aquifer was confirmed at the TCE contaminated site at Lewisburg, Tennessee indicating that natural attenuation should not be disregarded (Byl et

al., 2002). Also, Randrianarivelo et al. conducted remedial investigations of a karst aquifer in Pennsylvania. They selected chemical oxidation in conjunction with monitored natural attenuation or land use control as preferred remediation method alternative based on the hydrogeological conditions of the site and the extent of the delineated contamination plumes (Randrianarivelo et al., 2017).

Regarding ISCO remediation method, by assuming that the contaminated area is well identified and the injection fluid has the right dosage and residence time, the possibility of delivering injection fluid to the contaminated area with minimal error is the major challenge, similar to other fluid-based remediation methods in karst aquifers. Moreover, for achieving the highest efficiency in treating contaminants diffused into the rock matrix or moving with slow advective transport, oxidizing agents should remain in the contaminated area. However, rapid movement of water through preferential flow pathways dilutes the oxidizing agents. Thus, it can be asserted that application of ISCO in karst aquifers is limited. Regardless of limitations in applying treatment methods for karst aquifers, Randrianarivelo et al pointed out chemical oxidation as a beneficial techniques for GW remediation (Randrianarivelo et al., 2017).

The main challenge associated with these techniques is to identify the zone that requires treatment. Based on site conditions and the type of contaminants, the most effective technology should be employed. Assessment of remediation technique requires an appropriate monitoring approach (for locations consist of springs, streams, extraction systems, and previously tested wells) and comprehensive hydrogeological and water quality sampling data (Randrianarivelo et al., 2017).

3.2.2 Remediation by mitigating exposure pathways Exposure pathways often play a crucial role in spreading the contamination through karst aquifers. Steele-Valentín and Padilla delineated the potential exposure pathway of contaminants in Vega Alta, Puerto Rico. The pathway starts from superfund sites and is directed toward Atlantic Ocean. Moreover, it covers a noticeable lateral extension which impacts wells, wetlands and streams in the area (Steele-Valentín and Padilla, 2009). Because of interconnectivity between conduits and fractures within karst aquifer, controlling exposure pathways is more difficult compared to other aquifer types. However, mitigating exposure pathways can be done by treating at the tap, replacing drinking water supplies, treating spring flow using active and passive methods, land cover control using fences, signage, deed restriction and local law enforcement. In spite of their high treatment capabilities, these techniques often require long-term operation and maintenance costs (Randrianarivelo et al., 2017). One of the promising, cost-effective methods for minimizing the exposure pathways is electrochemical treatment: environmentally friendly (option of using solar power), independent (no other remedial actions needed) and controllable with regard to the rates of redox reactions. The high removal percentage for various contaminants in karst groundwater has been achieved in laboratory and pilot scale studies (Fallahpour, 2016; Gregor et al., 2017; Mao et al., 2012).

4 Application of Modeling in Karst In order to predict the behavior aquifers based on hydrological variations, groundwater models have been developed by hydrogeologists and water resources scientists. In addition, some models are developed to chemically analyze the water quality and to simulate fate and transport of contaminants. A groundwater flow model is able to exhibit precise

representation of hydrological and geological systems and also it can give a real insight into relationship and interactions between system elements. Modeling using computer programs can be truly beneficial when there are karst aquifers in the case study location. This is mainly due to the fact that karst aquifers are very heterogeneous and anisotropic and have a complex structure. Ergo, developing a relatively sophisticated model is the best option for simulating these types of aquifers. 4.1 Model Parameters and Development Depending on the soil type and water table level, the percolation rate regarding the movement of water from saturated zone to groundwater differs (Ritchey and Rumbaugh, 1996). Additionally, the impact of human interferences with natural water cycle which can be caused mostly by irrigation and pumping water from wells, should be taken into account. Actually, a groundwater model can determine how much it is perilous for an aquifer and also for the ecosystem if certain level of human interference exists. This can help developing water resources management plans that can not only help optimizing water extraction, but also can preserve the environment and natural resources (Drew and Hötzl, 1999). Other physical parameters such as topographical and geological information of the region that is going to be modeled should be given to the modeling platform (Peterson and Wicks, 2006). Anisotropy of aquifers regarding hydraulic conductivity, which is a parameter that can have a dissimilar value in each direction, can only be considered in two or three-dimensional models. Nowadays, by developing computer programs and also due to the need of acquiring more valid results as the output of modeling, three-dimensional models are more acceptable despite the possible complex procedure of setting them up (Anderson et al., 2015). In fact, while onedimensional models can be applied for vertical flow in multiple horizontal layers and two-dimensional models considers water flow in a vertical plain and this is repeated in multiple parallel vertical plains as well, a three-dimensional model subdivides the flow region into smaller cells with possible different properties regarding aquifer condition, soil characteristics and water flow (Ebrahim, 2014). Mostly, numerical analysis and tools should be used to solve complex differential equations of groundwater flow. In fact, a groundwater flow model is able to represent conceptual model of an aquifer mathematically and this mathematical representation enables researchers to solve the governing equations numerically by computers (Ebrahim, 2014; World Meteorological Organization, 2009). Using numerical solutions for solving groundwater flow equations in a threedimensional scale is beneficial for models which follow the flow domain discretization approach. Usually, in a groundwater flow model, hydraulic head at each cell center and groundwater flow rate between cells can be considered as outcomes. Moreover, impacts on streamflow because of pumping or long-term impacts of current pumping can be assessed. Checking the consistency of datasets and parameters and also defining e framework for future studies related to groundwater and aquifer condition are two prominent products of a groundwater model (Hartmann et al., 2014; Ritchey and Rumbaugh, 1996). Researchers and hydrogeologists have tried to develop groundwater models which can predict and simulate the groundwater flow in the karst aquifers in a regional or local scale. Regional groundwater models are capable of optimizing groundwater resources development plans, analyzing water budget of aquifers and assessing regional flow systems. Zhou and Li published a review paper regarding regional groundwater modeling and discussed their characteristics and associates drawbacks (Zhou and Li, 2011). Sauter assessed the quantification and prediction of

regional groundwater flow and transport in a karst aquifer in Germany. He discussed how the most appropriate modeling tool can be selected and how it can be used for simulation for a specific case study location. By analyzing spring flow, spatiotemporal variations of groundwater levels and hydraulic parameters, the model was able to successfully describe the karst aquifer system in the studied location (Sauter, 1992). Figure 6 depicts the schematic diagram of the process of developing a groundwater model.

Figure 6. Schematic diagram representing the process of developing a groundwater model – Modified from (World Meteorological Organization, 2009) 4.2 Spatially lumped models and distributed parameter models Ghasemizadeh et al. categorized groundwater models into different groups based on their capabilities and characteristics. Due to high level of heterogeneity and anisotropy in karst aquifers, accurately understanding of their behavior and distribution has always been challenging. This has forced modelers to employ approximate-based approaches and consequently consider the impact of the uncertainties caused by these approaches in their models. Hence, Spatially Lumped Models (SLM) and Distributed Models (DM) or Spatiotemporal Distributed Models (SDM) were introduced as two general approaches in modeling karst aquifers (Ghasemizadeh et al., 2012; World Meteorological Organization, 2009). Spatially lumped models comprises of concentrated elements at spatially singular points; whereas, the elements are spatially distributed in distributed models. Hence, in distributed systems, physical quantities are spatially and temporally dependent. Spatially lumped (or global) models do not consider spatial alternation of flow patterns and are supposed to simulate a global chemical-hydrological response at the aquifer output point (for example spring discharge point) with regard to inputs of the aquifer (e.g. rivers, groundwater recharge points, net runoff etc.) (Ghasemizadeh et al., 2012; Singh, 2014). Assessing temporal alternations is an approach that spatially lumped models take to describe the global water balance and hydrological behavior of an aquifer. Moreover, in spatially lumped models, some factors that cause

complexity in calculations and simulating are neglected due to simplifying assumptions and hence, using only the global parameters in simple ordinary linear differential equations and also low data requirements, are some of their properties that can be considered when trying to select the best modeling approach for groundwater flow and transport simulation. Although these models cannot produce accurate results, especially in karstic areas, they have been widely used by researchers in the areas where less data is available or only the prediction of groundwater flow, spring discharge and groundwater levels is necessary (Long, 2015b; Panagopoulos, 2012). Hydrograph-Chemograph Analysis (Dewandel et al., 2003), Linear Storage Models (or Rainfall-Discharge Models) (Butscher and Huggenberger, 2008) and Soft Computing Techniques such as Fuzzy Logic (Mohd Adnan et al., 2013; Rezaei et al., 2013), Genetic Algorithm (McKinney and Lin, 1994; Nicklow et al., 2010) and Artificial Neural Network (ANN) (Hu et al., 2008), are three main approaches with regard to spatially lumped models that have been adopted by hydrological modelers. In contrast, distributed models take complex parameters involved in groundwater flow and transport into account. In these models, dependent hydrological parameters and boundary conditions can be spatiotemporally variable and this will require the equations to be solved numerically and based on partial differential equations (Asher et al., 2015; Kuniansky, 2016). Also, due to the fact that all variables should be defied to the system, collecting more data and paying careful attention to details in this type of modeling is demanded which can make it more challenging. (Dong et al., 2012; Long and Gilcrease, 2009). For karst aquifer modeling, different type of distributed models based on the level of simplified assumptions have been used. Each of the developed models and methods treats complexity of karst aquifer differently and simulates groundwater flow based on its own logic and assumptions. Equivalent Porous Medium (EPM), Double Porosity (or Continuum) Method (DPM), Discrete Fracture Network (DFN), Discrete Channel (or Conduit) Network (DCN) and Hybrid Models (HM) are five common modeling approaches in distributed systems that have their own characteristics which will be discussed shortly in the following (Ghasemizadeh et al., 2012). DFN can be subcategorized into Discrete singular fracture set approach (DSFS) and Discrete multiple fracture set (DMFS) approaches. Sometimes, EPM and DPM are also mentioned as Single continuum porous equivalent approach (SCPE) and Double Continuum porous equivalent approach (DCPE) respectively in the literature. It is worth mentioning that employing Hybrid Models, which are the result of integrating discrete models and EPM approach and are also called coupled continuum pipe flow models, can be beneficial in many cases regarding modeling complex hydrological systems such as karst aquifers (Kiraly, 1998; Liedl et al., 2003). Figure 7 demonstrated schematic configuration of the aforementioned distributed modeling approaches.

Figure 7. Distributed parameter modeling methods for karst aquifers – Modified from (Kuniansky, 2016) Bauer et al. have developed a numerical model to describe the influence of exchange flow between conduits and fissured system. They found out that under conditions of early karst evolution, conduit development is faster. Hence, exchange flow plays an important role in developing early karst evolution in limestone aquifers (Bauer et al., 2003). Also, some researchers have employed numerical modeling approaches to describe groundwater flow and transport in rough fractures (Briggs et al., 2014) and karst aquifers (Faulkner et al., 2009). However, modeling karst aquifers cannot only be carried out by numerical approaches (Barrett and Charbeneau, 1998). Furthermore, for simulating the genesis of karst aquifer systems, a numerical couple reactive network model, comprising of a 2D porous continuum flow module, a discrete pipe network for modelling flow and transport in the conduits and a carbonate dissolution module was developed by (Clemens T. , 1997). 4.3 Computer Models and Programs MODFLOW is the most common groundwater modeling code that has been used due to its capability of simulating complex groundwater flows in a three-dimensional scale. Working based on finite difference method and block-centered approach, MODFLOW simulates the groundwater within the aquifer by considering different type of layers underground (i.e. confined, unconfined or both) and also different recharge or discharge sources such as areal recharge, groundwater flow to wells, runoff caused by rainfall, flow to riverbeds, spring flow etc. (Harbaugh, 2005). The initial version of MODFLOW (MODFLOW-2000) was released in the year 2000 and five years later, the updated version (MODFLOW2005) started gaining attentions from groundwater modelers and hydrogeologists. To enhance the application of MODFLOW-2000, two models have been introduced by USGS which are VSF and MF2K-GWT. Basically, VSF is a version of MODFLOW-2000 that in addition to the ability of MODFLOW-2000 to model groundwater flow using a finite-difference method in a 3-D scale, can be applicable for variably saturated flow (VSF) (Thoms et al., 2006). MF2K-GWT is an integrated model with MODFLOW-2000 that have the ability to simulate groundwater flow and solute transport (U.S. Geological Survey Website, 2012). Nevertheless, some programs that were independent to MODFLOW but developed by USGS were released as well such as HST3D (3-D Heat and Solute Transport Model) that is able to simulate ground-water flow and associated heat and solute transport in a 3D scale. Its capabilities can be used in analyzing problems associates with landfill leaching, seawater intrusion, hot-water geothermal systems etc. (Kipp, 1997). The most updated version of MODFLOW program (MODFLOW 6) was released recently. In this program, any number of models can be used for concurrent simulation. These models can have inter-connection with each other and this can help solving complex hydrogeological problems in many cases such as the conditions in karst aquifers. Also, within this framework, multiple local GW models can be coupled with regional scale models (Langevin et al., 2017). Moreover, Conduit Flow Package (CFP), which can be coupled with MODFLOW-2005, can facilitate simulation of karstic geometry and GW movement and consequently, increase the accuracy of GW flow modeling in conduits (Shoemaker et al., 2007). After releasing MODFLOW-2005, several associated models and packages were introduced and released based on numerous approaches and techniques. As an example, MT3D model, which is a modular, comprehensive, numerical three-dimensional solute transport model, was developed by USGS. This model has been designed to work very well

regarding simulation of solute transport and reactive solute transport in complex hydrological systems. Being connected to MODFLOW, which is the USGS groundwater flow simulator, MT3D is able to simulate and analyze advectiondominated transport, especially solute transport, without refining new models (Bedekar et al., 2016). Lautz and Siegel used MT3D and MODFLOW to simulate groundwater and surface water mixing in the hyporheic zone. They took advantage of this model due to its ability to simulate advective transport and source and sink mixing of solutes (Lautz and Siegel, 2006). Taking advantage of the features in MODFLOW and MT3D, a new computer program, SEAWAT, was released to assist hydrogeologists in simulating three-dimensional, variable-density and transient groundwater flow that can be coupled with solute transport. In the last version of SEAWAT (version 4), the effect of fluid viscosity and density fluctuations can be considered in simulation of groundwater flow and solute transport. This will allow the users to recognize this model as a tool that can be used in a wide range of simulation practices including seawater intrusion in coastal aquifers (Langevin, 2009; Langevin et al., 2008). Xu employed SEAWAT in his dissertation to study seawater intrusion into a coastal karstic aquifer in Florida and achieved accurate results. It is worth mentioning that seawater intrusion can be considered as a substantial source of brackish water in coastal aquifers such as karst aquifer in north coast of Puerto Rico (Xu, 2016). In addition, FEFLOW (Finite Element Subsurface Flow System) is a finite-element package for simulating 3D and 2D fluid density-coupled flow, contaminant mass (salinity) and heat transport in the subsurface. It has several applications including regional groundwater management, saltwater intrusion, seepage through dams and levees, land use and climate change scenarios, groundwater remediation and natural attenuation and also groundwater-surface water interaction. As an example, a study has been conducted to simulate groundwater dynamics in an irrigation and drainage network in Uzbekistan using FEFLOW. After model calibration and validation, the results show high level of accuracy and can be used for hydrogeological management plans (Diersch, 2014; Khalid Awan et al., 2015). SUTRA is another model that has been released for simulating 2-D saturated-unsaturated, fluid-density-dependent flow with energy transport or chemically-reactive single-species solute transport capable of analyzing saltwater intrusion and energy transport. It uses a 2D hybrid finite-element and integrated finite-difference approach to approximate the governing flow and transport equations that explain the two interdependent processes. It should be noted that the 3D version of this model has also been released recently. In SUTRA’s Version 2.2 specification of time-dependent boundary conditions can be identified without programming FORTRAN code. SUTRA, can also describe chemical species transport including absorption, production and decay processes and assess well performance and pumping test data (Voss and Provost, 2002). For instance, Hussain et al. used SUTRA in their paper to study coastal aquifer systems that are subjected to seawater intrusion (Hussain et al., 2015). Visual MODFLOW Flex platform, an integrated modeling environment that connects MODFLOW and MT3D, is able to simulate complex 3D groundwater flow and contaminant transport. Its graphical user interface and 3D visualization capabilities in addition to its ability to simulate groundwater flow and contaminant transport can gain attention of hydrological and groundwater modelers. As an example of work, Varghese, Raikar and Purandara successfully developed a Visual MODFLOW Flex model for simulation of groundwater flow in a region in India (Kumar and Singh, 2015; Varghese et al., 2015).

CHEMFLO-2000, which is interactive software for simulating water and chemical movement in unsaturated soils, enables users to simulate groundwater flow and chemical fate and transport in vadose zones. The model can be used as a tool that can enhance the understanding of unsaturated flow and transport processes. In this model, water movement and chemical transport are modeled using the Richards and the convection-dispersion equations, respectively. The equations are solved numerically using the finite differences approach (Nofziger and Wu, 2003). Another 3D finite-element based model for simulating flow and transport is 3DFEMFAT. This model works for saturated/unsaturated heterogeneous and anisotropic media. Its typical applications include infiltration, agriculture pesticides, sanitary landfill, hazardous waste disposal sites, density-induced flow and transport, saltwater intrusion, etc. Its flexibility and feasibility in simulating a wide range of practical problems especially by employing its transport module, has made it valuable software for researchers and transport modelers. Also its application in studying seawater intrusion in coastal aquifer has been verified by some scientists (Lathashri and Mahesha, 2016; Park et al., 2012). Regarding surface water and groundwater interaction which was discussed in the previous sections in detail, GSFLOW (Groundwater and Surface-water FLOW) was released by USGS in 2008 as an integrated tool that is able to couple groundwater and surface water flow models by taking advantage of the approaches used in USGS Precipitation-Runoff Modeling System (PRMS) and the USGS Modular Groundwater Flow Model (MODFLOW and MODFLOW-NWT). Meteorological and hydrological data such as rainfall, sunny hours and temperature in addition to groundwater stresses and initial/boundary conditions are involved as inputs for the process of simulation in this model. GSFLOW can also take into account the impact of land cover change, climate change and groundwater extraction on surface water and groundwater flow for spatiotemporally variable situations (Markstrom et al., 2008). However, regarding its limitations, it was asserted by researchers that its ability to simulate surface water and groundwater in karst aquifers with high level of heterogeneity is not guaranteed (Fulton et al., 2015). By taking advantage of a control volume finite-difference method, MODFLOW-USG (Un-Saturated Grid version of MODFLOW) is able to simulate groundwater flow and its related processes. This version of MODFLOW supports different types of structured and unstructured grids. This capability is extremely useful when high resolution along rivers and around wells is needed. In addition, MODFLOW-USG couples Connected Linear Network (CLN) process to Groundwater Flow (GWF) process, which was introduced in MODFLOW-2005, to analyze and simulate the influence of karst conduits and multi-node wells. Hence, this version can help modelers to gain a deeper understanding about karst systems and conduit networks (Panday et al., 2013). Moreover, for the purpose of generating layered quadtree grids that can be used in MODFLOW-USG or other similar numerical models, a new computer program, GRIDGEN, was developed by Lien et al. in 2015. After reading a 3-D base grid, GRIDGEN continues dividing into refinement features, which has already been provided by user, until reaching the desired refinement level. After finishing the process of gridding, a tree structure file will be created and can be used in numerical models such as MODFLOW-USG. This model was used for assessing the Biscayne aquifer in southern Florida in which karst aquifers are abundant (Lien et al., 2014). Developing a Newton-Raphson Formulation for MODFLOW-2005 for offering an enhanced solution for problems related to groundwater flow in unconfined aquifers, MODFLOW-NWT has been introduced and developed by (Niswonger et al., 2011). Its main application in addition to Surface-Water Routing (Hughes et al., 2012) and Seawater Intrusion (Bakker et al., 2013) can be described as its ability to solve problems that are coupled with drying and rewetting nonlinearities in equations that govern groundwater flow in unconfined aquifers.

A recently developed model similar to MODFLOW but with a wider range of applicability in describing hydrological systems is Rainfall-Response Aquifer and Watershed Flow Model (RRAWFLOW). This lumped-parameter model receives hydrological inputs such as rainfall, recharge and discharge etc. and is able to simulate groundwater level, streamflow and spring flow. It also can be used for modeling solute transport in aquifers and assessing system response to hydrological events (Long, 2015a). For classification of karst aquifers and characterizing time-variant systems, Long and Mahler developed and used this model in 2013. This model was used to predict and classify hydraulic responses to recharge in two karst aquifers in Texas and South Dakota, USA (Long and Mahler, 2013). Usually, groundwater flow and contaminant transport models are used simultaneously using software platforms such as GMS. Several researchers have done flow and transport analysis (e.g. using MODFLOW and MT3D) and have achieved accurate and valid results (Abdalla and Khalaf, 2015; Bora and Borah, 2016). Also, few scientists have studied the groundwater flow and contaminant transport in karstic aquifer of northern Puerto Rico using GMS and their modeling results show its capability in analyzing and describing hydrological systems with complex properties such as high level of heterogeneity and anisotropy (Ghasemizadeh, 2015; Ghasemizadeh et al., 2016; Maihemuti et al., 2015). Table 2 elaborates the characteristics and application of aforementioned most commonly used groundwater modeling codes. Table 2. Prevalent groundwater models that have been used for simulating groundwater flow and contaminant transport. FE and FD represent Finite Element and Finite Difference respectively.

Heat Transport

Solute Transport

GW Flow

Model/Software

Modeling technique

Focus

3DFEMFAT FE

*

*

AQUA3D FE

*

*

FD

*

*

*

CHEMFLO

FEFLOW FE

*

*

*

GSFLOW FD

*

HST3D FD

*

FD

*

MODFLOW

*

*

Application and Advantages / Reference

GW modeling in saturated/unsaturated heterogeneous and anisotropic media, simulation of infiltration, agriculture pesticides, sanitary landfill, hazardous waste disposal sites, density-induced flow and transport, seawater intrusion etc.(Lathashri and Mahesha, 2016; Park et al., 2012) 3D groundwater flow and transport simulation for homogeneous and anisotropic flow conditions, simulation of heat and contaminant transport by taking into account the effect of dispersion Simulation of water movement and chemical fate and transport in vadose zones and layered soil by employing improved numerical methods (Nofziger and Wu, 2003) regional groundwater management, saltwater intrusion, seepage through dams and levees, land use and climate change scenarios, groundwater remediation and natural attenuation, groundwater-surface water interaction (Diersch, 2014; Khalid Awan et al., 2015) Coupled Groundwater and surface water model which can assess the hydrological behavior based on land use change, climate variability and groundwater withdrawals (Markstrom et al., 2008) sub-surface-waste injection, landfill leaching, saltwater intrusion, freshwater recharge and recovery, radioactivewaste disposal, hot water geothermal systems, and subsurface-energy storage (Kipp, 1997) Simulation of steady or unsteady flow in complex flow system with irregular geometry, Simulation of flow from external stresses in a confined or unconfined aquifer

Ease of Use

Accuracy for Karst Modeling

High

Medium

Medium to High

Low

Low

Medium

High

Medium

Medium

Low

Medium

Medium

High

Medium to High

MODFLOWNWT

FD

*

MODFLOWOWHM FD

MODFLOWUSG

FD

*

*

MT3D FD

*

SEAWAT FD

*

*

*

SURTA FE

*

*

*

(Harbaugh, 2005), High applicability for karst aquifers if it couples with CFP package Surface water and groundwater interactions, seawater intrusion and solving problems related to drying and rewetting nonlinearities of the unconfined GW flow equation (Niswonger et al., 2011) Simulation, analysis, and management of human and natural water movement within a physically-based supplyand-demand framework, seawater intrusion, conjunctive use of groundwater and surface water (Hanson et al., 2014) Unstructured grid version of MODFLOW for simulating GW flow and other related processes, simulation of the effects of multi-node wells, karst conduits and tile drains (Panday et al., 2013) simulation of solute transport and reactive solute transport in complex hydrological systems and analyzing advection-dominated solute transport (Bedekar et al., 2016) 3D simulation of variable density, transient groundwater flow in porous media coupled with multi-species solute and heat transport, seawater intrusion in coastal aquifers (Langevin, 2009; Langevin et al., 2008; Post, 2011) Simulation of saturated-unsaturated, fluid-densitydependent groundwater flow with energy transport or chemically-reactive single-species solute transport (Voss and Provost, 2002)

Medium

Low

Low Low to Medium

High

Medium to High

High

Medium to High

High

Medium

Medium

Low to Medium

4.4 Equivalent Porous Media (EPM) method Several approaches have been followed to achieve accurate and valid results with acceptable efficiency at the same time. For some cases, simulation of groundwater hydraulic and contaminant transport in karst aquifers has been carried out by employing Equivalent Porous Media (EPM) method (Scanlon et al., 2003). Basically, using EPM approach for modeling a karst aquifer means considering simplifying assumptions in order to make the model more practical and applicable. Ghasemizadeh et al. developed their model based in EPM approach and found out that its result is acceptable for predicting water table fluctuations. Although their EPM-based model was not expected to be accurate enough for contaminant transport, they found good agreement between their model output and actual data regarding spreading TCE contaminant (Ghasemizadeh et al., 2015). Furthermore, in another study, by employing drainage features in regional groundwater flow modeling in karstic aquifer of northern Puerto Rico, Ghasemizadeh et al. asserted that they were able to improve their simulation by assigning arrays of adjacent model cells with drains to simulate conduits. They suggested that using this feature can be truly helpful especially when there is not sufficient data for conduit characteristics. Similarly, Maihemuti et al. developed a MODFLOW model using EPM approach for assessing karst aquifer system and groundwater resources for a case study location in northern Puerto Rico. Their model was supposed to predict the karst system response to rainfall events and high pumping demands and also to describe the hydrological behavior of the existing aquifer. They concluded that although there is high potential of conduit dominated flow, the result of their EPMbased approach is reliable in representing the hydrodynamics of the karst aquifer in their case study location (Maihemuti et al., 2015). 4.5 How Remote Sensing Can Improve Karst GW Modeling? Remote sensing tools, either coupled with modeling methods or be used separate, can be very beneficial because of their strong data analysis and management capability which allows assessment of several datasets and layers simultaneously.

However, lack of high resolution data for local studies can be problematic in some cases (Manda and Gross, 2006; Theilen-Willige et al., 2014). Using Geographic Information System (GIS) as a tool in groundwater modeling procedure, can be truly helpful; because all parameters such as distribution of rainfall, groundwater recharge and discharge and also land cover are defined within a spatial context (Singh and Fiorentino, 1996). Several researchers have took advantage of this powerful tool directly or as a parallel method in integrated approaches (Dar et al., 2010; Nampak et al., 2014). (Alonso-Contes, 2011) used remote sensing and advanced digital image processing techniques to delineate karst features which can enhance the understanding with regard to hydrogeology of the Tanamá River and Rio Grande de Arecibo catchments located in the north coast tertiary basin of Puerto Rico. Basically, remote sensing tools have assisted the author in lineament mapping for GW exploration. Also, Manda and Gross have employed GIS analysis to characterize solution conduits in karstic areas. Based on their study, they have shown that GIS-based methods can be used for determining depths, dimensions, shapes, apertures and connectivity of potential conduits and also for describing physical characteristics that have an effect on the groundwater flow in karst aquifers (Manda and Gross, 2006). In addition, Theilen-Willige et al. employed GIS and remote sensing methods by analyzing satellite data in order to detect of near-surface faults and fracture zones that can lead to dissolution processes in conduits of karst aquifers (Theilen-Willige et al., 2014). Numerous researchers have taken advantage of remote sensing techniques to improve their models and solve some unanswered and complex problems by increasing the accuracy of prediction and also by taking into account other hydrological phenomena (Ashraf and Ahmad, 2012; Machiwal et al., 2012; Thakur et al., 2016; Xu et al., 2011). Table 3 elaborates the potential application of GIS and remote sensing in different phases of integrated groundwater modeling. Table 3. The potential role of GIS and remote sensing in different steps of groundwater modeling procedure – Modified from (Ashraf and Ahmad, 2012) Phase

GIS functions

Modeling Steps

Data Collection and Analysis

Data input, Digitization, Data conversion (import/export, Coordinate transformation, Map retrieval

Groundwater and hydrological data collection

Developing Conceptual Model

Conversion of vector and raster layers, Data integration, Image processing, buffering, Surface generation, Linking of spatial and attribute data

Developing conceptual model

Model Design

Map calculations, Neighborhood operations, Interpolation, Theissen polygons, buffering, Surface generation

Delineating boundary conditions, Mesh generation, 3D layering of the aquifer

Data layers integration

Parameter zonation, Recharge estimation, Water balance

Overlay analysis

Steady-state and Transient-state simulations

Statistical analysis

Parameters estimation

Data retrieval

Prediction, Assessing different scenarios

Data visualization, Presentation of simulated results

Map composition

Model Calibration

Model Generalization, Predictions and Result Presentation

5 Conclusion Employing a comprehensive and efficient approach for managing water resources in regions where groundwater is the key source of water supply and vulnerable to contamination (e.g. karst aquifers) is vital. Limestone karstic aquifer of northern Puerto Rico has been experiencing high level of contamination and this has resulted in accelerating rate of preterm births in the island in addition to other human health related disorders. In PR, several wells and sites were considered as the National Priority List (NPL) Superfund Sites and remediation action for treating water in these sites are in process. By considering this region as the case study location and also by focusing on karst aquifers as the most complicated and hard-to-analyze forms of aquifers, a review study was presented to qualitatively and quantitatively assess GW resources. After an explanation of karst systems and their associated study methods in addition to surface water and GW interaction (SWGWI), a brief review discussion on groundwater contamination and risk assessment was presented. Potential contamination threats in karst aquifers were discussed and different remediation techniques were evaluated. Furthermore, a comprehensive discussion on existing groundwater modeling methods with regard to their application, advantages and disadvantages was presented. Mostly, numerical modeling approaches are taken by modelers to simulate complex groundwater systems. Numerical modeling tools often take advantage of finite-difference and finite-element techniques to solve complicated equations that governs hydrological system dynamics in an aquifer. For karst aquifers, lumped models and spatially distributed models are considered as two general modeling approaches that can be used for certain conditions. In addition, various computer-based models such as MODFLOW have been explained and evaluated. A combined set of models can be used to simulate groundwater flow and contaminant transport either in steady state or transient form. Additionally, application of remote sensing and GIS in assessing water resources in karst and also for modeling purposes was explained. Overall, it can be asserted that karst aquifers have the most complex behavior compared to other aquifer types. Because of their high level of permeability, they are always prone to contamination. Hence, employing different methods for studying karst aquifer and their inter-relationship with surface water and ground surface may not always lead to accurate results. However, although integrated techniques (taking advantage of different methods for a case study) and also collecting sufficient data can be costly, they can greatly enhance the understanding of a certain karst aquifer. These integrated techniques and data collection processes are strongly recommended for north coast karst aquifer of Puerto Rico where there are several superfund sites and human health and ecosystem are at high risk. From sustainable development perspective, environmental, economic and social impacts are consequences of water pollution in any region of the world. Hence, careful attention should be paid to preserve water resources. References Abdalla, M.G. , Khalaf, S., 2015. Groundwater Modeling of Multi-Aquifer Systems Using GMS. J. Wastewater. Treat. Analysis., 06 (01). doi:10.4172/2157-7587.1000184 Acosta-Colón, Á.A. 2016. Cave Characterization in the North Karst Belt Zone of Puerto Rico: Cave Mesofauna Diversity as an Indicator of Pathogenic and Opportunistic Species. In White, W.B.;Herman, J.S.;Herman, E.K. and Rutigliano, M. (Eds.), Karst Groundwater Contamination and Public Health: Springer. Alloway, B.J. 2013. Sources of Heavy Metals and Metalloids in Soils Heavy Metals in Soils (Third Edition ed., Vol. 22). Dordrecht: Springer. Almasri, M.N., 2007. Nitrate contamination of groundwater: A conceptual management framework. Environ. Impact. Assess. Rev., 27 (3), 220-242. doi:10.1016/j.eiar.2006.11.002 Almasri, M.N. , Kaluarachchi, J.J., 2004. Assessment and management of long-term nitrate pollution of ground water in agriculturedominated watersheds. J. Hydrol., 295 (1-4), 225-245. doi:10.1016/j.jhydrol.2004.03.013

Alonso-Contes, C. 2011. Lineament mapping for groundwater exploration using remotely sensed imagery in a karst terrain: Rio Tanama and Rio de Arecibo basins in the northern karst of Puerto Rico. (Masters), Michigan Technological University. Anderson, M.P., Woessner, W.W. , Hunt, R.J., 2015. Applied Groundwater Modeling: Simulation of Flow and Advective Transport (Second Edition) Elsevier Inc. Angulo, B., Morales, T., Uriarte, J.A. , Antigüedad, I., 2011. Hydraulic conductivity characterization of a karst recharge area using water injection tests and electrical resistivity logging. Eng. Geol., 117 (1-2), 90-96. doi:10.1016/j.enggeo.2010.10.008 Arfib, B., de Marsily, G. , Ganoulis, J., 2007. Locating the zone of saline intrusion in a coastal karst aquifer using springflow data. Ground Water, 45 (1), 28-35. doi:10.1111/j.1745-6584.2006.00252.x Asher, M.J., Croke, B.F.W., Jakeman, A.J. , Peeters, L.J.M., 2015. A review of surrogate models and their application to groundwater modeling. Water Resour. Res., 51 (8), 5957-5973. doi:10.1002/2015wr016967 Ashraf, A. , Ahmad, Z. 2012. Integration of Groundwater Flow Modeling and GIS. In Nayak, P. (Ed.), Water Resources Management and Modeling (pp. 239-262): InTech. ASTM,1995. D 5717–95: Guide for Design of Ground-Water Monitoring Systems in Karst and Fractured Rock Aquifers, ASTM International Babiker, I., 2004. Assessment of groundwater contamination by nitrate leaching from intensive vegetable cultivation using geographical information system. Environ. Int., 29 (8), 1009-1017. doi:10.1016/s0160-4120(03)00095-3 Bailly-Comte, V., Jourde, H. , Pistre, S., 2009. Conceptualization and classification of groundwater–surface water hydrodynamic interactions in karst watersheds: Case of the karst watershed of the Coulazou River (Southern France). 376 (3-4), 456-462. doi:10.1016/j.jhydrol.2009.07.053 Bakalowicz, M., 2005. Karst groundwater: a challenge for new resources. Hydrogeol. J., 13 (1), 148-160. doi:10.1007/s10040-0040402-9 Bakker, M., Schaars, F., Hughes, J.D., Langevin, C.D. , Dausman, A.M. 2013. Documentation of the Seawater Intrusion (SWI2) Package for MODFLOW, U.S. Geological Survey, Techniques and Methods 6-A46 Barrett, M.E. , Charbeneau, R.J., 1998. A parsimonious model for simulating flow in a karst aquifer. J. Hydrol. , 196 (1-4), 47-65. doi:doi.org/10.1016/S0022-1694(96)03339-2 Baskaran, S., Ransley, T., Brodie, R.S. , Baker, P., 2009. Investigating groundwater–river interactions using environmental tracers. Aust. J. Earth Sci., 56 (1), 13-19. doi:10.1080/08120090802541887 Bauer-Gottwein, P., Gondwe, B.R.N., Charvet, G., Marín, L.E., Rebolledo-Vieyra, M. , Merediz-Alonso, G., 2011. Review: The Yucatán Peninsula karst aquifer, Mexico. Hydrogeol. J., 19 (3), 507-524. doi:10.1007/s10040-010-0699-5 Bauer, S., Liedl, R. , Sauter, M., 2003. Modeling of karst aquifer genesis: Influence of exchange flow. Water Resour. Res., 39 (10), 112. doi:10.1029/2003wr002218 Bayless, R., Cinotto, P.J., Ulery, R.L., Taylor, C.J., McCombs, G.K., Kim, M.H. , Nelson, H.L. 2014. Surface-water and karst groundwater interactions and streamflow-response simulations of the karst-influenced upper Lost River watershed, Orange County, Indiana, Scientific Investigations Report Bedekar, V., Morway, E.D., Langevin, C.D. , Tonkin, M. 2016. MT3D-USGS version 1: A U.S. Geological Survey release of MT3DMS updated with new and expanded transport capabilities for use with MODFLOW, U.S. Geological Survey, Techniques and Methods 6-A53 Beyke, G.L., Hodges, B.A. , Jones, G.N., 2014. Electrical Resistance Heating of Volatile Organic Compounds in Sedimentary Rock. Remediation J., 25 (1), 53-70. doi:10.1002/rem.21414 Biaggi, N., 1995. Puerto Rico's Water Pollution Image. Water Pollut. Control Federation, 37 (3), 381-391. Bonneau, J., Fletcher, T.D., Costelloe, J.F. , Burns, M.J., 2017. Stormwater infiltration and the ‘urban karst’ – A review. J. Hydrol., 552, 141-150. doi:10.1016/j.jhydrol.2017.06.043 Bora, G. , Borah, T., 2016. Simulation of Flow and Transport Processes of a Non-uniform Aquifer by GMS. J. Civ. Eng. Environ. Tech., 3 (3), 237-240. Brian Katz, Coplen, T., Bullen, T. , Davis, H., 1997. Use of Chemical and Isotopic Tracers to Characterize the Interactions Between Ground Water and Surface Water in Mantled Karst. 35 (6), 1014-1028. doi:10.1111/j.1745-6584.1997.tb00174.x Briggs, S., Karney, B.W. , Sleep, B.E., 2014. Numerical modelling of flow and transport in rough fractures. J. Rock Mech. Geotech. Eng., 6 (6), 535-545. doi:10.1016/j.jrmge.2014.10.004 Brodie, R., Sundaram, B., Tottenham, R., Hostetler, S. , Ransley, T. 2007. An Overview of Tools for Assessing Groundwater-Surface Water Connectivity Dept. of Agriculture, Fisheries and Forestry, Australian Governemnt, Butscher, C. , Huggenberger, P., 2008. Intrinsic vulnerability assessment in karst areas: A numerical modeling approach. Water Resour. Res., 44 (3), 1-15. doi:10.1029/2007wr006277 Byl, T.D., Hileman, G.E., Williams, S.D. , Farmer, J.J. 2001. Geochemical and microbial evidence of fuel biodegradation in a contaminated karst aquifer in southern Kentucky, U.S. Geological Survey Karst Interest Group Proceedings, WaterResources Investigations Report 01-4011 Byl, T.D., Hileman, G.E., Williams, S.D., Metge, D.W. , Harvey, R.W. 2002. Microbial Strategies for Degradation of Organic Contaminants in Karst. Paper presented at the U.S. Geological Survey Artificial Recharge Workshop, Sacramento, California. Byl, T.D. , Williams, S.D. 2000. Biodegradation of Chorinated ethenes at a karst site in Middle Tennessee, U.S. Geological Survey, Water-Resources Investigations Report 99-4285 Carmona De Jesús, M. , Padilla., I.Y. 2015. Characterization of TCE NAPL and Dissolved Phase Transport in Karst Media. Paper presented at the American Geophysical Union 2015 Fall Meeting, San Francisco, CA. Castro-Prieto, J., Martinuzzi, S., Radeloff, V.C., Helmers, D.P., Quiñones, M. , Gould, W.A., 2017. Declining human population but increasing residential development around protected areas in Puerto Rico. Biol. Conserv., 209, 473-481. doi:10.1016/j.biocon.2017.02.037

Chalikakis, K., Plagnes, V., Guerin, R., Valois, R. , Bosch, F.P., 2011. Contribution of geophysical methods to karst-system exploration: an overview. Hydrogeol. J., 19 (6), 1169-1180. doi:10.1007/s10040-011-0746-x Chen, H., Liu, J., Zhang, W. , Wang, K., 2011. Soil hydraulic properties on the steep karst hillslopes in northwest Guangxi, China. Environ. Earth Sci., 66 (1), 371-379. doi:10.1007/s12665-011-1246-y Chen, Z., Auler, A.S., Bakalowicz, M., Drew, D., Griger, F., Hartmann, J., Jiang, G., Moosdorf, N., Richts, A., Stevanovic, Z., Veni, G. , Goldscheider, N., 2017. The World Karst Aquifer Mapping project: concept, mapping procedure and map of Europe. Hydrogeol. J., 25 (3), 771-785. doi:10.1007/s10040-016-1519-3 Cherry, G.S. 2001. Simulation of Flow in the Upper North Coast Limestone Aquifer, Manatí-Vega Baja Area, Puerto Rico, U.S. Geological Survey, Water-Resources Investigations Report 00-4266 Chu, H., Wei, J., Wang, R. , Xin, B., 2016. Characterizing the interaction of groundwater and surface water in the karst aquifer of Fangshan, Beijing (China). 25 (2), 575-588. doi:10.1007/s10040-016-1507-7 Clemens T. , H.D., Sauter M. , Liedl R. , Teutsch G. 1997. Modelling the genesis of karst aquifer systems using a coupled reactive network model. Paper presented at the Hard Rock Geosciences (Proceedings of Rabat Symposium S2). Conboy, M.J. , Goss, M.J., 2000. Natural protection of groundwater against bacteria of fecal origin. J. Contam. Hydrol, 43 (1), 1-24. Conde-Costas, C. , Gómez-Gómez, F. 1999. Assessment of nitrate contamination of the upper aquifer in the Manati-Vega Baja area, Puerto Rico, U.S. Geological Survey, Water-Resources Investigations Report 99-4040 Dar, I.A., Sankar, K. , Dar, M.A., 2010. Remote sensing technology and geographic information system modeling: An integrated approach towards the mapping of groundwater potential zones in Hardrock terrain, Mamundiyar basin. J. Hydrol., 394 (3-4), 285-295. doi:10.1016/j.jhydrol.2010.08.022 Deeb, R., Hawley, E., Kell, L. , O’Laskey, R. 2011. Assessing Alternative Endpoints for Groundwater Remediation at Contaminated Sites, Environmental Security Technology Certification Program (ESTCP), ESTCP Project Report ER-200832 Dewandel, B., Lachassagne, P., Bakalowicz, M., Weng, P. , Al-Malki, A., 2003. Evaluation of aquifer thickness by analysing recession hydrographs. Application to the Oman ophiolite hard-rock aquifer. J. Hydrol., 274 (1-4), 248-269. doi:10.1016/s0022-1694(02)00418-3 Diersch, H.-J.G., 2014. FEFLOW finite element modeling of Flow, Mass and Heat Transport in Porous and Fractured Media: Springer. Dodgen, L.K., Kelly, W.R., Panno, S.V., Taylor, S.J., Armstrong, D.L., Wiles, K.N., Zhang, Y. , Zheng, W., 2017. Characterizing pharmaceutical, personal care product, and hormone contamination in a karst aquifer of southwestern Illinois, USA, using water quality and stream flow parameters. Sci. Total Environ., 578, 281-289. doi:10.1016/j.scitotenv.2016.10.103 Dong, Y., Li, G. , Xu, H., 2012. An areal recharge and discharge simulating method for MODFLOW. Comput. Geosci., 42, 203-205. doi:10.1016/j.cageo.2011.10.005 Drew, D. , Hötzl, H., 1999. Karst Hydrogeology and Human Activities: Impacts, Consequences and Implications: CRC Press. Ebrahim, G.Y., 2014. Modelling Groundwater Systems: Understanding and Improving Groundwater Quantity and Quality Management: CRC Press. El Alfy, M. , Faraj, T., 2017. Spatial distribution and health risk assessment for groundwater contamination from intensive pesticide use in arid areas. Environ. Geochem. Health, 39 (1), 231-253. doi:10.1007/s10653-016-9825-1 Engelhardt, I., De Aguinaga, J.G., Mikat, H., Schuth, C. , Liedl, R., 2014. Complexity vs. simplicity: groundwater model ranking using information criteria. Ground Water, 52 (4), 573-583. doi:10.1111/gwat.12080 Fallahpour, N. 2016. Effect of natural organic matter, metal ions, and nitrate on electrochemical dechlorination of trichloroethylene. (PhD Dissertation), Northeastern University. Faulkner, J., Hu, B.X., Kish, S. , Hua, F., 2009. Laboratory analog and numerical study of groundwater flow and solute transport in a karst aquifer with conduit and matrix domains. J. Contam. Hydrol., 110 (1-2), 34-44. doi:10.1016/j.jconhyd.2009.08.004 Fetter, C.W., 2001. Applied HydroGeology 4th Edition: Prentice Hall, Inc. Field, M.S. 2017. Investigating and Remediating Contaminated Karst Aquifers Karst Groundwater Contamination and Public Health (pp. 101-115): Springer. Fleckenstein, J.H., Krause, S., Hannah, D.M. , Boano, F., 2010. Groundwater-surface water interactions: New methods and models to improve understanding of processes and dynamics. Adv. Water Resour., 33 (11), 1291-1295. doi:10.1016/j.advwatres.2010.09.011 Fleury, P., Bakalowicz, M. , de Marsily, G., 2007. Submarine springs and coastal karst aquifers: A review. J. Hydrol., 339 (1-2), 7992. doi:10.1016/j.jhydrol.2007.03.009 Ford, D. , Williams, P.D., 2007. Karst Hydrogeology and Geomorphology: John Wiley & Sons Ltd. Freeze, A. , Cherry, J., 1979. Groundwater: Prentice Hall, Inc. Fu, T., Chen, H., Zhang, W., Nie, Y. , Wang, K., 2015. Vertical distribution of soil saturated hydraulic conductivity and its influencing factors in a small karst catchment in Southwest China. Environ. Monit. Assess., 187 (3), 92. doi:10.1007/s10661015-4320-1 Fulton, J.W., Risser, D.W., Regan, R.S., Walker, J.F., Hunt, R.J., Niswonger, R.G., Hoffman, S.A. , Markstrom, S. 2015. WaterBudget and Recharge-Area Simulations for Spring Creek and Nittany Creek Basins and Parts of the Spruce Creek Basin, Centre and Huntingdon Counties, Pennsylvania, Water Years 2000–06, U.S. Geological Survey, Scientific Investigations Report 2015-5073 Gárfias-Soliz, J., Llanos-Acebo, H. , Martel, R., 2009. Time series and stochastic analyses to study the hydrodynamic characteristics of karstic aquifers. Hydrol. Process., 24, 300-316. doi:10.1002/hyp.7487 Ghasemizadeh, R. 2015. Modeling Groundwater flow and Contaminant transport in the North Coast Limestone karst aquifer system of Puerto Rico. (PhD Dissertation), Northeastern University, Boston, MA. Ghasemizadeh, R., Hellweger, F., Butscher, C., Padilla, I., Vesper, D., Field, M. , Alshawabkeh, A., 2012. Review: Groundwater flow and transport modeling of karst aquifers, with particular reference to the North Coast Limestone aquifer system of Puerto Rico. Hydrogeol. J., 20 (8), 1441-1461. doi:10.1007/s10040-012-0897-4

Ghasemizadeh, R., Yu, X., Butscher, C., Hellweger, F., Padilla, I. , Alshawabkeh, A., 2015. Equivalent Porous Media (EPM) Simulation of Groundwater Hydraulics and Contaminant Transport in Karst Aquifers. PLoS One, 10 (9), 1-21. doi:10.1371/journal.pone.0138954 Ghasemizadeh, R., Yu, X., Butscher, C., Padilla, I.Y. , Alshawabkeh, A., 2016. Improved regional groundwater flow modeling using drainage features: a case study of the central northern karst aquifer system of Puerto Rico (USA). Hydrogeol. J., 24 (6), 1463-1478. doi:10.1007/s10040-016-1419-6 Giudici, M., Margiotta, S., Mazzone, F., Negri, S. , Vassena, C., 2012. Modelling hydrostratigraphy and groundwater flow of a fractured and karst aquifer in a Mediterranean basin (Salento peninsula, southeastern Italy). Environ. Earth Sci., 67 (7), 1891-1907. doi:10.1007/s12665-012-1631-1 Gleeson, T., Moosdorf, N., Hartmann, J. , Beek, L.P.H.v., 2014. A glimpse beneath earth's surface: GLobal HYdrogeology MaPS (GLHYMPS) of permeability and porosity. Geophys. Res. Lett., 41 (11). doi:10.1002/2014GL059856 Goldscheider, N. , Drew, D., 2007. Methods in Karst Hydrogeology. London: Taylor & Francis. Goldscheider, N., Meiman, J., Pronk, M. , Smart, C., 2008. Tracer tests in karst hydrogeology and speleology. Int. J. Speleol., 37 (1), 27-40. Gómez-Gómez, F., Rodríguez-Martínez, J. , Santiago, M. 2014. Hydrogeology of Puerto Rico and the outlying islands of Vieques, Culebra, and Mona, U.S. Geological Survey, Scientific Investigations Map 3296 González-Pinzón, R., Ward, A.S., Hatch, C.E., Wlostowski, A.N., Singha, K., Gooseff, M.N., Haggerty, R., Harvey, J.W., Cirpka, O.A. , Brock, J.T., 2015. A field comparison of multiple techniques to quantify groundwater–surface-water interactions. Freshw. Sci., 34 (1), 139-160. doi:10.1086/679738 Gregor, S., Fallahpour, N., Rajic, L. , Alshawabkeh, A. 2017. Electrochemical Remediation of Contaminated Groundwater: Pilot Scale Study. In White, W.;Herman, J.;Herman, E. and Rutigliano, M. (Eds.), Karst Groundwater Contamination and Public Health (pp. 117-120): Springer. Guay, C., Nastev, M., Paniconi, C. , Sulis, M., 2013. Comparison of two modeling approaches for groundwater-surface water interactions. Hydrol. Process., 27 (16), 2258-2270. doi:10.1002/hyp.9323 Hanson, R.T., Boyce, S.E., Schmid, W., Hughes, J.D., Mehl, S.M., Leake, S.A., III, T.M. , Niswonger, R.G. 2014. One-Water Hydrologic Flow Model (MODFLOW-OWHM), U.S. Geological Survey, Techniques and Methods 6-A51 Harbaugh, A.W. 2005. MODFLOW-2005, The U.S. Geological Survey Modular Ground-Water Model—the Ground-Water Flow Process, U.S. Geological Survey Techniques and Methods 6-A16 Hartmann, A., Goldscheider, N., Wagener, T., Lange, J. , Weiler, M., 2014. Karst water resources in a changing world: Review of hydrological modeling approaches. Rev. Geophys., 52 (3), 218-242. doi:10.1002/2013rg000443 Hashim, M.A., Mukhopadhyay, S., Sahu, J.N. , Sengupta, B., 2011. Remediation technologies for heavy metal contaminated groundwater. J. Environ. Manage., 92 (10), 2355-2388. doi:10.1016/j.jenvman.2011.06.009 Hennessey, M., Fischer, M. , Staples, J.E., 2016. Zika virus spreads to new areas—region of the Americas, May 2015–January 2016. Am. J. Transplant., 16 (3), 1031-1034. Hodges, B., Jones, G., Beyke, G. , Crownover, C. (2014). Thermal Remediation of Karst Limestone – Redstone Arsenal, Alabama Retrieved from http://www.thermalrs.com/news-story/59 Hu, B.X., Meerschaert, M.M., Barrash, W., Hyndman, D.W., He, C., Li, X. , Guo, L., 2009. Examining the influence of heterogeneous porosity fields on conservative solute transport. J. Contam. Hydrol., 108 (3-4), 77-88. doi:10.1016/j.jconhyd.2009.06.001 Hu, C., Hao, Y., Yeh, T.C.J., Pang, B. , Wu, Z., 2008. Simulation of spring flows from a karst aquifer with an artificial neural network. Hydrol. Process., 22 (5), 596-604. doi:10.1002/hyp.6625 Hughes, J.D., Langevin, C.D., Chartier, K.L. , White, J.T. 2012. Documentation of the Surface-Water Routing (SWR1) Process for Modeling Surface-Water Flow with the U.S. Geological Survey Modular Groundwater Model (MODFLOW–2005), U.S. Geological Survey, Techniques and Methods 4-A40 Hussain, M.S., Javadi, A.A., Ahangar-Asr, A. , Farmani, R., 2015. A surrogate model for simulation–optimization of aquifer systems subjected to seawater intrusion. J. Hydrol., 523, 542-554. doi:10.1016/j.jhydrol.2015.01.079 Illman, W.A., Liu, X. , Craig, A., 2007. Steady-state hydraulic tomography in a laboratory aquifer with deterministic heterogeneity: Multi-method and multiscale validation of hydraulic conductivity tomograms. J. Hydrol., 341 (3-4), 222-234. doi:10.1016/j.jhydrol.2007.05.011 Jankowski, J. 2007. Surface Water-Groundwater Interactions in a Catchment Impacted by Longwall Mining. Paper presented at the Triennial Conference on Mine Subsidence, Wollongong, Australia. Jones, I.C. , Banner, J.L., 2003. Estimating recharge thresholds in tropical karst island aquifers: Barbados, Puerto Rico and Guam. J. Hydrol., 278 (1-4), 131-143. doi:10.1016/s0022-1694(03)00138-0 Kačaroğlu, F., 1999. Review of Groundwater Pollution and Protection in Karst Areas. Water Air Soil Pollut., 113 (1), 337-356. doi:10.1023/A:1005014532330 Kalbus, E., Reinstorf, F. , Schirmer, M., 2006. Measuring methods for groundwater-surface water interactions: a review. Hydrol. Earth Syst. Sci., 10 (6), 873-887. Katz, B.G., Coplen, T.B., Bullen, T.D. , Davis, J.H., 1997. Use of Chemical and Isotopic Tracers to Characterize the Interactions Between Ground Water and Surface Water in Mantled Karst. Ground Water, 35 (6), 1014-1028. Kaufmann, G., 2003. Modelling unsaturated flow in an evolving karst aquifer. J. Hydrol., 276 (1-4), 53-70. doi:10.1016/s00221694(03)00037-4 Kaufmann, G., 2016. Modelling karst aquifer evolution in fractured, porous rocks. J. Hydrol., 543, 796-807. doi:10.1016/j.jhydrol.2016.10.049 Kendall, C. , Aravena, R. 2000. Nitrate Isotopes in Groundwater Systems. In Cook, P.G. and Herczeg, A.L. (Eds.), Environmental Tracers in Subsurface Hydrology (pp. 261-297).

Khalid Awan, U., Tischbein, B. , Martius, C., 2015. FEFLOW-Simulating Groundwater Dynamics Using Feflow-3D Groundwater Model Under Complex Irrigation and Drainage Network of Dryland Ecosystems of Central Asia. Irrig. Drain., 64, 283-296. doi:10.1002/ird.1897 Kipp, K.L. 1997. Guide to the Revised Heat and Solute Transport Simulator: HST3D - Version 2, U.S. Geological Survey, WaterResources Investigations Report 97-4157 Kiraly, L., 1998. Modelling karst aquifers by the combined discrete channel and continuum approach. Bull. du Centre d'hydrogéologie, 16, 77-98. Kiraly, L., 2003. Karstification and Groundwater Flow. Speleogenesis Evol. Karst Aquifers, 1 (3), 1-26. Kovács, A., Perrochet, P., Király, L. , Jeannin, P.Y., 2005. A quantitative method for the characterisation of karst aquifers based on spring hydrograph analysis. J. Hydrol., 303 (1-4), 152-164. doi:10.1016/j.jhydrol.2004.08.023 Kumar, C.P. , Singh, S., 2015. Concepts and Modeling of Groundwater System. Int. J. Innov. Res. Sci. Eng. Technol., 2 (2), 262-271. Kuniansky, E.L. 2016. Simulating Groundwater Flow in Karst Aquifers with Distributed Parameter Models— Comparison of PorousEquivalent Media and Hybrid Flow Approaches, Scientific Investigations Report 2016-5116 Langevin, C.D. 2009. SEAWAT: a computer program for simulation of variable-density groundwater flow and multi-species solute and heat transport, U.S. Geological Survey, Fact Sheet FS 2009-3047 Langevin, C.D., Daniel T. Thorne, J., Dausman, A.M., Sukop, M.C. , Guo, W. 2008. SEAWAT Version 4: A Computer Program for Simulation of Multi-Species Solute and Heat Transport, U.S. Geological Survey Techniques and Methods 6-A22 Langevin, C.D., Hughes, J.D., Banta, E.R., Niswonger, R.G., Panday, S. , Provost, A.M. 2017. Documentation for the MODFLOW 6 Groundwater Flow Model, Techniques and Methods Lapworth, D.J., Baran, N., Stuart, M.E. , Ward, R.S., 2012. Emerging organic contaminants in groundwater: A review of sources, fate and occurrence. Environ. Pollut., 163, 287-303. doi:10.1016/j.envpol.2011.12.034 Lathashri, U.A. , Mahesha, A., 2016. Predictive Simulation of Seawater Intrusion in a Tropical Coastal Aquifer. J. Environ. Eng., 142 (12), 1-13. doi:10.1061/(asce)ee.1943-7870.0001037 Lautz, L.K. , Siegel, D.I., 2006. Modeling surface and ground water mixing in the hyporheic zone using MODFLOW and MT3D. Adv. Water Resour., 29 (11), 1618-1633. doi:10.1016/j.advwatres.2005.12.003 Liedl, R., Sauter, M., Hückinghaus, D., Clemens, T. , Teutsch, G., 2003. Simulation of the development of karst aquifers using a coupled continuum pipe flow model. Water Resour. Res., 39 (3), 1-11. doi:10.1029/2001wr001206 Lien, J.M., Liu, G. , Langevin, C.D. 2014. GRIDGEN Version 1.0: A Computer Program for Generating Unstructured Finite-Volume Grids, U.S. Geological Survey, Open-File Report 2014–1109 Loáiciga, H.A. 2001. Ground-Water/Surface-Water Interactions in a Karst Aquifer. Paper presented at the Specialty Symposium on Integrated Surface and Ground Water Management at the World Water and Environmental Resources Congress, Orlando, Florida. https://ascelibrary.org/doi/abs/10.1061/40562%28267%2916 Long, A.J., 2015a. RRAWFLOW: Rainfall-Response Aquifer and Watershed Flow Model (v1.15). 8, 865-880. doi:10.5194/gmd-8865-2015 Long, A.J., 2015b. RRAWFLOW: Rainfall-Response Aquifer and Watershed Flow Model (v1.15). Geosci. Model Dev., 8 (3), 865880. doi:10.5194/gmd-8-865-2015 Long, A.J. , Gilcrease, P.C., 2009. A one-dimensional heat-transport model for conduit flow in karst aquifers. J. Hydrol., 378 (3-4), 230-239. doi:10.1016/j.jhydrol.2009.09.024 Long, A.J. , Mahler, B.J., 2013. Prediction, time variance, and classification of hydraulic response to recharge in two karst aquifers. Hydrol. Earth Syst. Sci., 17 (1), 281-294. doi:10.5194/hess-17-281-2013 Luo, M., Chen, Z., Criss, R.E., Zhou, H., Huang, H., Han, Z. , Shi, T., 2016. Dynamics and anthropogenic impacts of multiple karst flow systems in a mountainous area of South China. Hydrogeol. J., 24 (8), 1993-2002. doi:10.1007/s10040-016-1462-3 Machiwal, D., Mishra, A., Jha, M.K., Sharma, A. , Sisodia, S.S., 2012. Modeling Short-Term Spatial and Temporal Variability of Groundwater Level Using Geostatistics and GIS. Nat. Resour. Res., 21 (1), 117-136. doi:10.1007/s11053-011-9167-8 Mahler, B.J., Personné, J.C., Lods, G.F. , Drogue, C., 2000. Transport of free and particulate-associated bacteria in karst. J. Hydrol., 238 (3-4), 179-193. Maihemuti, B., Ghasemizadeh, R., Yu, X., Padilla, I. , Alshawabkeh, A.N., 2015. Simulation of Regional Karst Aquifer System and Assessment of Groundwater Resources in Manati-Vega Baja, Puerto Rico. J. Water Resour. Prot., 07 (12), 909-922. doi:10.4236/jwarp.2015.712075 Manda, A.K. , Gross, M.R., 2006. Identifying and characterizing solution conduits in karst aquifers through geospatial (GIS) analysis of porosity from borehole imagery: An example from the Biscayne aquifer, South Florida (USA). Adv. Water Resour., 29 (3), 383-396. doi:10.1016/j.advwatres.2005.05.013 Mao, X., Yuan, S., Fallahpour, N., Ciblak, A., Howard, J., Padilla, I., Loch-Caruso, R. , Alshawabkeh, A.N., 2012. Electrochemically induced dual reactive barriers for transformation of TCE and mixture of contaminants in groundwater. Environ. Sci. Technol., 46 (21), 12003-12011. doi:10.1021/es301711a March of Dimes Website. (2016). 2016 Premature Birth Report Cards. from March of Dimes http://www.marchofdimes.org/mission/prematurity-reportcard.aspx?utm_source=homeslide&utm_medium=MODsite&utm_campaign=reportcard16&utm_content=report-card# Markstrom, S.L., Niswonger, R.G., Regan, R.S., Prudic, D.E. , Barlow, P.M. 2008. GSFLOW—Coupled Ground-Water and SurfaceWater Flow Model Based on the Integration of the Precipitation-Runoff Modeling System (PRMS) and the Modular GroundWater Flow Model (MODFLOW-2005), U.S. Geological Survey, Techniques and Methods 6-D1 Martinez, J.L., Raiber, M. , Cox, M.E., 2015. Assessment of groundwater-surface water interaction using long-term hydrochemical data and isotope hydrology: Headwaters of the Condamine River, Southeast Queensland, Australia. Sci. Total Environ., 536, 499-516. doi:10.1016/j.scitotenv.2015.07.031 Mathes, M.V. 2000. Relation of bacteria in limestone aquifers to septic systems in Berkeley County, West Virginia, U.S. Geological Survey, Water-Resources Investigations Report 2000-4229

Mathews, T. , MacDorman, M. 2011. Infant mortality statistics from the 2007 period linked birth/infant death data set, Centers for Disease Control and Prevention (CDC), National Vital Statistics Report 59(6) McKinney, D.C. , Lin, M.-D., 1994. Genetic algorithm solution of groundwater management models. Water Resour. Res., 30 (6), 1897-1906. Metcalfe, C.D., Beddows, P.A., Bouchot, G.G., Metcalfe, T.L., Li, H. , Van Lavieren, H., 2011. Contaminants in the coastal karst aquifer system along the Caribbean coast of the Yucatan Peninsula, Mexico. Environ. Pollut., 159 (4), 991-997. doi:10.1016/j.envpol.2010.11.031 Meyerhoff, S.B., Maxwell, R.M., Revil, A., Martin, J.B., Karaoulis, M. , Graham, W.D., 2013. Characterization of groundwater and surface water mixing in a semiconfined karst aquifer using timelapse electrical resistivity tomography. Water Resour. Res., 50, 2566-2585. doi:10.1002/2013WR013991 Milanović, S., Dragišić, V., Radulović, M.M. , Stevanović, Z. 2015. Prevent Leakage and Mixture of Karst Groundwater Karst Aquifers—Characterization and Engineering (pp. 531-599): Springer. Mohd Adnan, M.R.H., Sarkheyli, A., Mohd Zain, A. , Haron, H., 2013. Fuzzy logic for modeling machining process: a review. Artif. Intell. Rev., 43 (3), 345-379. doi:10.1007/s10462-012-9381-8 Mongelli, G., Monni, S., Oggiano, G., Paternoster, M. , Sinisi, R., 2013. Tracing groundwater salinization processes in coastal aquifers: a hydrogeochemical and isotopic approach in the Na-Cl brackish waters of northwestern Sardinia, Italy. Hydrol. Earth Syst. Sci., 17 (7), 2917-2928. doi:10.5194/hess-17-2917-2013 Morales, T., Angulo, B., Uriarte, J.A., Olazar, M., Arandes, J.M. , Antiguedad, I., 2017. Solute transport characterization in karst aquifers by tracer injection tests for a sustainable water resource management. J. Hydrol., 547, 269-279. doi:10.1016/j.jhydrol.2017.02.009 Morasch, B., 2013. Occurrence and dynamics of micropollutants in a karst aquifer. Environ. Pollut., 173, 133-137. doi:10.1016/j.envpol.2012.10.014 Muldoon, M.A., Borchardt, M.A., Spencer, S.K., Hunt, R.J. , Owens, D. 2016. Using Enteric Pathogens to Assess Sources of Fecal Contamination in the Silurian Dolomite Aquifer: Preliminary Results. In White, W.B.;Herman, J.S.;Herman, E.K. and Rutigliano, M. (Eds.), Karst Groundwater Contamination and Public Health: Springer. Nampak, H., Pradhan, B. , Manap, M.A., 2014. Application of GIS based data driven evidential belief function model to predict groundwater potential zonation. J. Hydrol., 513, 283-300. doi:10.1016/j.jhydrol.2014.02.053 Nicklow, J., Reed, P., Savic, D., Dessalegne, T., Harrell, L., Chan-Hilton, A., Karamouz, M., Minsker, B., Ostfeld, A., Singh, A. , Zechman, E., 2010. State of the Art for Genetic Algorithms and Beyond in Water Resources Planning and Management. Water Resour. Plan. Manage., 136 (4), 412-432. doi:10.1061/共ASCE兲WR.1943-5452.0000053 Niswonger, R.G., Panday, S. , Ibaraki, M. 2011. MODFLOW-NWT, A Newton Formulation for MODFLOW-2005, U.S. Geological Survey, Techniques and Methods 6–A37 Nofziger, D.L. , Wu, J. 2003. User's Manual: CHEMFLO-2000, Interactive Software for Simulating Water and Chemical Movement in Unsaturated Soils U.S. Environmental Protection Agency, Report EPA/600/R-03/008 Padilla, I., Irizarry, C. , Steele., K., 2011. Historical contamination of groundwater resources in the north coast karst aquifer of Puerto Rico. Rev. Dimens., 3 (7-12). Panagiotakis, I. , Dermatas, D., 2017. Soil and Groundwater Contamination and Remediation. Bull. Environ. Contam. Toxicol., 98 (3), 297-298. doi:10.1007/s00128-017-2043-4 Panagopoulos, G., 2012. Application of MODFLOW for simulating groundwater flow in the Trifilia karst aquifer, Greece. Environ. Earth Sci., 67 (7), 1877-1889. doi:10.1007/s12665-012-1630-2 Panday, S., Langevin, C.D., Niswonger, R.G., Ibaraki, M. , Hughes, J.D. 2013. MODFLOW–USG Version 1: An Unstructured Grid Version of MODFLOW for Simulating Groundwater Flow and Tightly Coupled Processes Using a Control Volume FiniteDifference Formulation, U.S. Geological Survey, Techniques and Methods 6-A45 Parise, M., Ravbar, N., Živanović, V., Mikszewski, A., Kresic, N., Mádl-Szőnyi, J. , Kukurić, N. 2015. Hazards in Karst and Managing Water Resources Quality. In Stevanović, Z. (Ed.), Karst Aquifers—Characterization and Engineering (pp. 601687): Springer. Park, S.U., Kim, J.M., Yum, B.W. , Yeh, G.T., 2012. Three-Dimensional Numerical Simulation of Saltwater Extraction Schemes to Mitigate Seawater Intrusion due to Groundwater Pumping in a Coastal Aquifer System. J. Hydrol. Eng., 17 (1), 10-22. doi:10.1061/(asce)he.1943-5584.0000412 Peterson, E.W. , Wicks, C.M., 2006. Assessing the importance of conduit geometry and physical parameters in karst systems using the storm water management model (SWMM). J. Hydrol., 329 (1-2), 294-305. doi:10.1016/j.jhydrol.2006.02.017 Post, V.E.A., 2011. A new package for simulating periodic boundary conditions in MODFLOW and SEAWAT. Comput. Geosci., 37 (11), 1843-1849. doi:10.1016/j.cageo.2011.01.012 Quinn, J.J., Tomasko, D. , Kuiper, J.A., 2006. Modeling complex flow in a karst aquifer. Sediment. Geol., 184 (3-4), 343-351. doi:10.1016/j.sedgeo.2005.11.009 Rafael, M.T., Ronald, R. , Anastacio, E., 2016. Analysis of Groundwater in Puerto Rico. Am. J. Water Resour., 4 (3), 68-76. doi:10.12691/ajwr-4-3-3 Randrianarivelo, M., Zhou, W. , Barsa, M., 2017. Remedial investigations of karst aquifers: a case study at former Marietta Air Force Station, Lancaster County, Pennsylvania. Carbonates Evaporite. doi:10.1007/s13146-017-0369-y Rezaei, F., Safavi, H.R. , Ahmadi, A., 2013. Groundwater vulnerability assessment using fuzzy logic: a case study in the Zayandehrood aquifers, Iran. Environ. Manage., 51 (1), 267-277. doi:10.1007/s00267-012-9960-0 Ritchey, J. , Rumbaugh, J., 1996. Subsurface Fluid-Flow (Ground-Water and Vadose Zone) Modeling: ASTM International. Rivera, V.L., Padilla, I.Y. , Torres, N.I. 2017. Spatiotemporal Assessment of CVOC Contamination in Karst Groundwater Sources and Exposure at Tap Water Point of Use. In White, W.;Herman, J.;Herman, E. and Rutigliano, M. (Eds.), Karst Groundwater Contamination and Public Health (pp. 231-235): Springer.

Rodriguez-Martinez, J. 1995. Hydrogeology of the North Coast Limestone Aquifer System of Puerto Rico, Water-Resources Investigations Report 94-4249 Rosenberry, D.O. , LaBaugh, J.W. 2008. Field Techniques for Estimating Water Fluxes Between Surface Water and Ground Water, U.S. Geological Survey Techniques and Methods 4-D2 Rugel, K., Golladay, S.W., Jackson, C.R. , Rasmussen, T.C., 2016. Delineating groundwater/surface water interaction in a karst watershed: Lower Flint River Basin, southwestern Georgia, USA. 5, 1-19. doi:10.1016/j.ejrh.2015.11.011 Rutigliano, M. 2016. Beyond Case Reports: Placing Karst in Context in Public Health Response to Groundwater Contamination. Paper presented at the Karst Groundwater Contamination and Public Health Conference, San Juan, Puerto Rico Sauter, M. 1992. Quantification and forecasting of regional groundwater flow and transport in a karst aquifer (Gallusquelle, Malm, SW. Germany). (PhD Dissertation), Universität Tübingen. Retrieved from http://hdl.handle.net/10900/48838 Scanlon, B.R., Mace, R.E., Barrett, M.E. , Smith, B., 2003. Can we simulate regional groundwater flow in a karst system using equivalent porous media models? Case study, Barton Springs Edwards aquifer, USA. J. Hydrol., 276 (1-4), 137-158. doi:10.1016/s0022-1694(03)00064-7 Shiklomanov, I.A. 1993. World fresh water resources. In Gleick, P.H. (Ed.), Water in Crisis. New York. Shoemaker, B., Kuniansky, E.L., Birk, S., Bauer, S. , Swain, E.D. 2007. Documentation of a Conduit Flow Process (CFP) for MODFLOW-2005, Techniques and Methods, Singh, A., 2014. Groundwater resources management through the applications of simulation modeling: a review. Sci. Total Environ., 499, 414-423. doi:10.1016/j.scitotenv.2014.05.048 Singh, V. , Fiorentino, M., 1996. Geographical Information Systems in Hydrology. Netherlands: Springer Sinreich, M. 2014. Contaminant Attenuation in Karst Aquifers: A Paradigm Shift. In Mudry, J.;Zwahlen, F.;Bertrand, C. and LaMoreaux, J. (Eds.), H2Karst Research in Limestone Hydrogeology (pp. 175-184): Springer. Smith, L.A., Barbour, S.L., Hendry, M.J., Novakowski, K. , van der Kamp, G., 2016. A multiscale approach to determine hydraulic conductivity in thick claystone aquitards using field, laboratory, and numerical modeling methods. Water Resour. Res., 52 (7), 5265-5284. doi:10.1002/2015wr018448 Sreekanth, J. , Datta, B., 2015. Review: Simulation-optimization models for the management and monitoring of coastal aquifers. Hydrogeol. J., 23 (6), 1155-1166. doi:10.1007/s10040-015-1272-z Steele-Valentín, K.M. , Padilla, I.Y. 2009. Assessment of Potential Exposure Pathways in Karst Groundwater Systems in Vega Alta, Puerto Rico using Geographic Information Systems. Paper presented at the 7th LACCEI Latin American and Caribbean Conference for Engineering and Technology, San Cristóbal, Venezuela. https://www.scribd.com/document/294991216/Assessment-of-Potential-Exposure-Pathways-in-Karst-Vega-Alta-PR Stevanović, Z., 2015. Karst Aquifers – Characterization and Engineering: Springer. Sudicky, E.A., Illman, W.A., Goltz, I.K., Adams, J.J. , McLaren, R.G., 2010. Heterogeneity in hydraulic conductivity and its role on the macroscale transport of a solute plume: From measurements to a practical application of stochastic flow and transport theory. Water Resour. Res., 46 (1), n/a-n/a. doi:10.1029/2008wr007558 Sui, Q., Cao, X., Lu, S., Zhao, W., Qiu, Z. , Yu, G., 2015. Occurrence, sources and fate of pharmaceuticals and personal care products in the groundwater: A review. Emerg. Contam., 1 (1), 14-24. doi:10.1016/j.emcon.2015.07.001 Sutton, J.E., Screaton, E.J. , Martin, J.B., 2014. Insights on surface-water/groundwater exchange in the upper Floridan aquifer, northcentral Florida (USA), from streamflow data and numerical modeling. Hydrogeol. J., 23 (2), 305-317. doi:10.1007/s10040014-1213-2 Thakur, J.K., Singh, S.K. , Ekanthalu, V.S., 2016. Integrating remote sensing, geographic information systems and global positioning system techniques with hydrological modeling. Appl. Water Sci. doi:10.1007/s13201-016-0384-5 Theilen-Willige, B., Malek, H., Charif, A., El Bchari, F. , Chaïbi, M., 2014. Remote Sensing and GIS Contribution to the Investigation of Karst Landscapes in NW-Morocco. Geosci., 4 (2), 50-72. doi:10.3390/geosciences4020050 Thoms, R.B., Johnson, R.L. , Healy, R.W. 2006. User's guide to the Variably Saturated Flow (VSF) Process for MODFLOW, U.S. Geological Survey, Techniques and Methods 6-A18 Thorn, R.H. , Coxon, C.E., 1992. Hydrogeological Aspects of Bacterial Contamination Some Western Ireland Karstic Limestone Aquifers. Environ. Geol. Water Sci., 20 (1), 65-72. U.S. Geological Survey Website. (2012). MF2K-GWT: Three-dimensional ground-water flow and solute-transport model integrated with MODFLOW-2000. Retrieved from https://water.usgs.gov/nrp/gwsoftware/mf2k_gwt/mf2k_gwt.html USEPA. 2003. Environmental Indicator Groundwater Contamination Chevron Phillips Chemical Puerto Rico Core LLC, Guayama PR, Environmental Protection Agency, Documantation of Environmental Indicator Determination USEPA. 2012. Summary of Technical Impracticability Waivers at National Priorities List Sites United States Environmental Protection Agency, OSWER Directive Report 9230.2-24 USEPA. 2013. Superfund Remedy Report (14th Edition), United Stated Environmental Protection Agency, Report EPA 542-R-13-016 Varghese, A.A., Raikar, R.V. , Purandara, B.K., 2015. Simulation of Groundwater Levels in Malaprabha Command Area using Visual MODFLOW Flex. Int. Res. J. Eng. Tech., 2 (8), 434-440. Vesper, D.J., Loop, C.M. , White, W.B., 2003. Contaminant transport in karst aquifers. Speleogenesis Evol. Karst Aquifers, 2 (1), 111. Vidal Montes, R., Martinez-Graña, A.M., Martínez Catalán, J.R., Ayarza Arribas, P. , Sánchez San Román, F.J., 2016. Vulnerability to groundwater contamination, SW salamanca, Spain. J. Maps, 12 (sup1), 147-155. doi:10.1080/17445647.2016.1172271 Voss, C.I. , Provost, A.M. 2002. SUTRA, A model for saturated-unsaturated variable-density ground-water flow with solute or energy transport, U.S. Geological Survey, Water-Resources Investigations Report 02-4231 Wakida, F.T. , Lerner, D.N., 2005. Non-agricultural sources of groundwater nitrate: a review and case study. Water Res., 39 (1), 3-16. doi:10.1016/j.watres.2004.07.026 Weary, D.J. , Doctor, D.H. 2014. Karst in the United States: a digital map compilation and database, U.S. Geological Survey, OpenFile Report 2014–1156

White, W.B. 2016. Contaminant Transport in Karst Aquifers: Systematics and Mechanisms. In White, W.B.;Herman, J.S.;Herman, E.K. and Rutigliano, M. (Eds.), Karst Groundwater Contamination and Public Health: Springer. Winter, T.C., Harvey, J.W., Franke, O.L. , Alley, W.M. 1998. Ground Water And Surface Water: A Single Resource, U.S. Geological Survey, Circular Report 1139 World Meteorological Organization. 2009. Guide to Hydrological Practices Volume II: Management of Water Resources and Application of Hydrological Practices (6th ed.): World Meteorological Organization (WMO). Wu, B., Zheng, Y., Tian, Y., Wu, X., Yao, Y., Han, F., Liu, J. , Zheng, C., 2014. Systematic assessment of the uncertainty in integrated surface water-groundwater modeling based on the probabilistic collocation method. Water Resour. Res., 50 (7), 5848-5865. doi:10.1002/2014wr015366 Xu, X., Huang, G., Qu, Z. , Pereira, L.S., 2011. Using MODFLOW and GIS to Assess Changes in Groundwater Dynamics in Response to Water Saving Measures in Irrigation Districts of the Upper Yellow River Basin. Water Resour. Manage., 25 (8), 2035-2059. doi:10.1007/s11269-011-9793-2 Xu, Z. 2016. Data analysis and numerical modeling of seawater intrusion through conduit network in a coastal karst aquifer. (PhD Dissertation), Florida State University, Tallahassee, Florida. Xu, Z., Bassett, S.W., Hu, B. , Dyer, S.B., 2016. Long distance seawater intrusion through a karst conduit network in the Woodville Karst Plain, Florida. Sci. Rep., 6 (32235), 1-10. doi:10.1038/srep32235 Yao, Z., Li, J., Xie, H. , Yu, C., 2012. Review on Remediation Technologies of Soil Contaminated by Heavy Metals. Proced. Environ. Sci., 16, 722-729. doi:10.1016/j.proenv.2012.10.099 Yu, X., Ghasemizadeh, R., Padilla, I., Irizarry, C., Kaeli, D. , Alshawabkeh, A., 2015. Spatiotemporal changes of CVOC concentrations in karst aquifers: analysis of three decades of data from Puerto Rico. Sci. Total Environ., 511, 1-10. doi:10.1016/j.scitotenv.2014.12.031 Yu, X., Ghasemizadeh, R., Padilla, I.Y., Kaeli, D. , Alshawabkeh, A., 2016. Patterns of temporal scaling of groundwater level fluctuation. J. Hydrol., 536, 485-495. doi:10.1016/j.jhydrol.2016.03.018 Zack, A., Rodriguez-Alonso, T. , Roman-Mas, A. 1987. Puerto Rico Ground-Water Quality, U.S. Geological Survey, Open-File Report 87-0749 Zemann, M., Wolf, L., Grimmeisen, F., Tiehm, A., Klinger, J., Hotzl, H. , Goldscheider, N., 2015. Tracking changing X-ray contrast media application to an urban-influenced karst aquifer in the Wadi Shueib, Jordan. Environ. Pollut., 198, 133-143. doi:10.1016/j.envpol.2014.11.033 Zhou, Y. , Li, W., 2011. A review of regional groundwater flow modeling. Geosci. Front., 2 (2), 205-214. doi:10.1016/j.gsf.2011.03.003

Highlights -

Karst aquifers are extremely permeable and vulnerable to contamination High rate of pre-term births in Puerto Rico due to GW contamination Assessing different types of contaminants in Karst Source & Plume Remediation and Remediation by mitigating exposure pathways in Karst GW flow and contaminant fate and transport simulation methods