To protect karst spring water resources, catchments must be known. We have developed a method for correlating spring hydrographs with newly available, high‐resolution, satellite‐based Global Precipitation Measurement data to rapidly and remotely locate recharge areas. We verify the method using a synthetic comparison of ground‐based rain gage data with the satellite precipitation data set. Application to karst springs is proven by correlating satellite data with hydrographs from well‐known springs with published catchments in Europe and North America. Application to an unknown‐catchment spring in Pennsylvania suggests distant recharge, requiring a flow path that crosses topographic divides, as well as multiple lithologies, physiographic provinces, and tectonic boundaries. Although surprising, this latter result is consistent with published geologic/geophysical, monitoring well, and stream gage data. We conclude that the method has considerable potential to improve the speed and accuracy of catchment identification and hydrodynamic characterization, with applications to water resource protection and groundwater exploration, among others.

We have developed a method for determining the recharge area for any spring using satellite‐based precipitation data. These data are compared to changes in flow of the spring. That is, when there are surges in the spring flow, we look at satellite data to see where it has recently rained. This will have important applications in protecting spring water resources worldwide.

1 Introduction Groundwater from karst aquifers supplies nearly 25% of human population [Ford and Williams, 2007]. Most large springs on Earth are karstic, with many vital to cities, including Vienna, Rome, and San Antonio [Kresic and Stevanovic, 2010]. Spring protection in the face of population growth and climate change requires specific strategies [Hartmann et al., 2014]. However, karst aquifers are often geologically complex and hydraulically connected over long distances [Bakalowicz, 2005]. Therefore, delineation of spring catchments requires three‐dimensional geological analysis, spring hydrograph monitoring and water balances, hydrochemical/isotopic methods, and artificial tracer tests [Goldscheider and Drew, 2007]. Groundwater flow paths in karst (confirmed by dye‐tracing) sometimes exceed 100 km, illustrating the expansiveness of these groundwater resources and the challenges of delineation [Ford and Williams, 1989]. This work is motivated by a karst spring in Pennsylvania, USA, locally called “the Bubble.” This spring, plus several smaller and connected ones nearby, discharges 5 to 7 times the expected infiltration in its topographic watershed [Becher, 1991]. The spring hydrograph (Figure 1) displays discharge surges, while temperature, conductivity, and turbidity remain constant. Water temperature varies seasonally by 0.3°C but is 6 months out of phase with air temperature. This suggests pressure‐pulse surges which we correlate, using a novel method, with precipitation recorded by NASA's Global Precipitation Measurement (GPM) satellite network [NASA, 2016]. Our goal was to constrain the extent of the unknown Bubble recharge zone. This procedure is tested using springs with published catchments. These tests, and the surprising results from the Bubble, support the potential worldwide utility of this method for identifying recharge areas and their hydrodynamics. Figure 1 Open in figure viewer PowerPoint Comparison of precipitation data for a TRG (teal), GPM (purple), and the Bubble hydrograph (black). The dashed blue lines highlight discharge events with insignificant local precipitation. The blue arrows show local precipitation with no change in Bubble discharge.

2 Methods The GPM constellation measures precipitation using microwaves. One data set is the “Integrated Multi‐satellitE Retrievals for GPM Level 3 Final Run” [Huffman et al., 2015]. This provides precipitation maps between 60°N and 60°S with 0.1° × 0.1° spatial resolution and 30 min temporal resolution. Data collection began in March 2014, and the data are available with a 4 month lag for calibration against ground‐based precipitation gages. The GPM Final Run is considered the most accurate of any remotely sensed precipitation data set [Huffman et al., 2015]. To check the accuracy of the GPM data, we recorded precipitation for 8 months using a tipping rain gage (TRG) near the Bubble. Figure 1 overlays TRG data on the GPM record for the Bubble pixel. While the magnitudes of events differ (the gage consistently records greater values—probably due to the inherent averaging within pixels in the GPM data), the times of events match closely (correlation coefficient = 0.66 at p = 0.05). The Bubble hydrograph is shown on Figure 1. In general, surge times correlate poorly with local TRG data; some surges occur without local precipitation and vice versa (Figure 1). One explanation is a component of remotely applied head from precipitation/recharge or associated changes in atmospheric pressure. Remote sources are consistent with discharge greatly exceeding infiltration in the topographic watershed. From the GPM maps, we construct a precipitation time series for every pixel within the search area, which is compared to the spring hydrograph. We developed an algorithm for making this comparison (termed ECHO for Empirically Constrained Hydrologic Operation) as summarized in Figure 2. Figure 2 Open in figure viewer PowerPoint Flowchart for ECHO analysis. The grey steps are automated. ECHO starts with a hydrograph as discharge D versus time [D, t]. Multiple GPM precipitation (P) maps [X, Y, P], spanning the time range of the hydrograph, are downloaded from NASA. ECHO combines individual map files into a 4‐D block of [X, Y, t, P]. The hydrograph is subjected to event recognition, which identifies an event with index (i) based on fitting a regression line to a moving window, and flagging the last point beyond which inclusion of additional points produces a positive change in slope greater than Δ, and which persists for more than n subsequent windows, where Δ and n are user‐selectable. The event list [i, t] is iteratively compared to the [t, P] series for each pixel using time lags spanning a specified range. For each pixel, the time lag that produces the maximum number of matches between spring discharge and precipitation events (“hits”) is recorded, along with the associated number of hits for that pixel. These map data [X, Y, hits] are exported and displayed as contours or classed postings. The method does not yet comprehend snowmelt signals, so winter months are excised from the hydrograph. The prototype algorithm was applied to the Bubble and blind‐tested at three karst springs of differing magnitudes (as defined by Meinzer [1923]) at diverse geographic locations and climates. These springs are described in Table 1 and shown on Figure 3. Table 1. Springs Used to Test the ECHO Prototype Name Location Low‐Flow Discharge (m3/s) Average Discharge (m3/s) High‐Flow Discharge (m3/s) Magnitude (Based on Average) Catchment (km2) The Bubble [Becher and Root, 1981 this study] Boiling Springs, Pennsylvania, USA 0.5 0.7 1 2nd Unknown Sägebach [Chen and Goldscheider, 2014 Goldscheider, 2005 Kleinwalsertal, Vorarlberg, Austria 0.15 0.35 3.5 2nd 28 Areuse [Kiraly, 2003 Canton de Neuchâtel, Switzerland 4 28 42 1st 127 Big Spring [Imes et al., 2007 Vineyard and Feder, 1982 Van Buren, Missouri, USA 7 12 39 1st 1104 Figure 3 Open in figure viewer PowerPoint Chen et al., 2017 Locations of test springs (Table 1 ) on the World Karst Aquifer Map by

3 Results ECHO results for the Bubble are depicted in Figure 4a. On this and subsequent maps, the hit count for each pixel is normalized by the maximum hit count for all pixels and time lags across the tested region so contour maps (from triangulation with linear interpolation) can use a consistent color scale. For the Bubble, we analyzed the full GPM global coverage (4.32 million pixels). Under 200 pixels registered >14 hits, and all (0.0043% of the total) falls within the map limits of Figure 4a. Two pixels registered the peak hit count of 21 (a ratio of approximately one in two million). The peaks occur in adjoining pixels approximately 60 km SE of the spring. The implied hydrogeological connection is surprising because the identified recharge area lies on the opposite side of the quartzite/metavolcanic ridge or “South Mountain anticlinorium” [Cloos, 1950] shown in Figure 4b. Figure 4 Open in figure viewer PowerPoint Stoeser et al., 2005 (a) ECHO output for the Bubble. The small black box is the Bubble GPM pixel. The colors depict hit count as a percentage of the peak which is 21 hits in two pixels east of W. Conewago. (b) Geology [from] and physiographic provinces (dashed). (c) Verification test showing GPM‐TRG hits at the Bubble. To confirm the ECHO result, Figure 4c shows a test using TRG data (Figure 1) as a mock hydrograph, correlated with the same GPM data as the Bubble test in Figure 4a. The highest hit count is in the Bubble/TRG pixel, validating the ECHO algorithm. Diminishing correlations outside the Bubble/TRG pixel are likely a semiquantitative indication of the resolution of the method at this location for this time range. That is, the lateral extent of precipitation events (beyond the Bubble pixel) and the typical storm track (SW to NE along the Great Valley) yield a spatial resolution depicted by the color contours of Figure 4c. We could use this pattern to correct for artificial enlargement of the recharge zone, but in the future, intend to prescreen the GPM data to retain only small fast‐moving storms. Figure 4c shows an important test of ECHO, but since the catchment for the Bubble is not known, this result cannot prove the effectiveness of the procedure. Therefore, ECHO was blind‐tested on three springs (Table 1 and Figure 3) with mapped catchments from traditional geologic mapping, water budget calculations, hydraulic modeling, and dye tracing confirmation. These include Sägebach (Austria), Areuse (Switzerland), and Big Spring (MO, USA). Figure 5 shows the geographic location for each spring, the published outline of the springshed, and the ECHO hit contours overlain on NASA Terra Moderate Resolution Imaging Spectroradiometer imagery. The ECHO algorithm used GPM data covering an area much larger than that in the figures, Eastern North America for Big Spring and all of Europe for Sägebach and Areuse. Figure 5 Open in figure viewer PowerPoint Contoured ECHO hits for (a) Areuse, (b) Sägebach, and (c) Big Spring. Catchments are blue, and springs are cyan. To characterize the effectiveness of ECHO and evaluate the significance of the contours in Figure 5, we performed a receiver operating characteristic (ROC) analysis [Hajian‐Tilaki, 2013] for the TGR test in Figure 4c and each of the blind tests in Figure 5 (see the supporting information). For each spring/test, the area enclosed by each successive hit count contour (from values of 100 through 92) was compared to the mapped catchment to calculate “detection completeness,” and “false alarm ratio.” Detection completeness is the fraction of the mapped catchment that a given hit count contour encloses (equivalent to specificity or true positive rate in general ROC analysis). False alarm ratio is the fraction of area enclosed by a contour that falls outside of the mapped catchment (equivalent to sensitivity or false alarm rate). The discrimination statistic d or area under the ROC curve [Nam and D'Agostino, 2002] was calculated for each test/spring and can be compared to theoretical “perfect” test and “worthless” test in the following list: perfect test, d = 1.00; TRG test, d = 0.91; Big Spring, d = 0.84; Sägebach, d = 0.73; Areuse, d = 0.63; and worthless test, d = 0.50. The TRG and spring applications show that ECHO clearly discriminates catchments at a statistically significant level. One limitation of this test is that it measures only the effectiveness of the method at matching the outline of the mapped catchment and does not value the ability of ECHO to simply locate each spring, with no additional input, in continent‐scale search areas. In addition, these discrimination statistics presume that the catchment maps are correct and that ECHO‐identified potential recharge zones outside the mapped catchment are actually “false positive.”

4 Discussion In each of the known‐catchment tests of ECHO, the greatest number of hits falls within the mapped recharge area, with counts falling‐off outside of it (similar to the GPM versus TRG test in Figure 4c) and good to excellent ROC behavior. An analysis of the time lag between precipitation and discharge peaks has not been completed and seems nonconstant for Areuse and Big Spring. For Sägebach, the consistent time lag matches dye trace transit times [Goldscheider, 2005; Göppert and Goldscheider, 2008]. For all springs, there are relatively high‐correlation pixels outside of the known springsheds (treated as “false positive” in the ROC analysis). This may be due to temporal coincidence of precipitation in these pixels with that in the springshed (unrelated to groundwater flow), or in some cases, ECHO might be identifying true additional, previously unidentified recharge areas. Longer time series comparisons with culled GPM data preserving only fast‐moving small‐extent storms could reduce spurious hits. Strong hits outside currently mapped springsheds (e.g., west of the marked Big Spring catchment) could motivate research to confirm surprising hydrogeological connections which often occur in karst terranes. Since there is no springshed map for the Bubble, and the ECHO‐located recharge area lies in a geologically unexpected location (Figure 4b), we compared the Bubble hydrograph to water data throughout mid‐Atlantic USA. We identified apparent correlation for two monitoring points shown on Figure 4a: Well AD‐146 [U.S. Geological Survey (USGS), 2016a] and W. Conewago Creek gaging station [USGS, 2016b]. Both display water‐level peaks that match surges in the Bubble (Figure 6). This across‐the‐mountain similarity of hydraulic behavior went unnoticed until the ECHO analysis. If ECHO is correct, this suggests that South Mountain may not be a hydrogeological barrier. An artificial tracer test could evaluate this hypothesis. A hydraulic fast‐path for the pressure pulses in the Bubble could be present if South Mountain is a crystalline allochthon that has overridden younger carbonates, consistent with the interpretation of seismic reflection data by Evans [1989] for the same “anticlinorium” farther south. A continuous karstified pathway, from the remote recharge area across the Gettysburg Basin and South Mountain, to the Bubble is possible since there are freshwater carbonates in the fault‐bounded rift basin [de Wet et al., 1998]. Detailed mapping and seismic studies are underway/planned to address these possibilities. We are also exploring the role of highly hydraulically conductive pseudo‐karstic fracture networks that may carry water through the noncarbonate rocks [Moser et al., 2014]. Figure 6 Open in figure viewer PowerPoint USGS, 2016a USGS, 2016b (top) Bubble hydrograph. (bottom) Piezometric level, Well AD‐146 [] and streamflow, W. Conewago Creek gaging station [] for a 2 week period in summer 2015. Well and gage locations shown in Figure 4 a. The current ECHO algorithm may not comprehend emergent discharge events with long lag times. This is especially true for springs that show large seasonal fluctuations. In addition, for existing spring discharge data—e.g., USGS National Water Information Mapper [USGS, 2016a]—the sampling intervals vary from hourly to daily or longer. ECHO handles this by comparing hydrograph event onset times to the 30 min sampled data. Resampling of the GPM and/or hydrograph data to a common interval will facilitate more rigorous correlation methods and should affect correlation behavior. This and other improvements to the methodology are underway.

5 Conclusions and Outlook Correlation of satellite GPM data with the Bubble hydrograph supports our hypothesis of a remote recharge component. Testing of the ECHO method using TRG data and hydrographs from other karst springs indicates that this will be a powerful tool with global applications for rapidly connecting water supplies to their recharge areas, whether local or remote (and unexpected). ECHO cannot replace hydrogeologic mapping, geochemical, geophysical, and dye tracing studies, but it could direct these to critical testing locations and provide a valuable groundwater exploration tool. In initial tests, ECHO has revealed a previously unrecognized correlation, and potential hydraulic connection, between water on either side of South Mountain in Pennsylvania. In addition, there are hints of unrecognized connections to Big Spring, Missouri. We anticipate that testing of other karst springs could reveal numerous new connections. We are currently working on improved methods for cross‐correlating satellite precipitation and hydrograph data to better‐constrain recharge locations and time lags and quantify associated error and significance statistics. We will include topographic filtering of cells with elevations below that of the discharge point. With these improvements, we intend to identify catchment areas and transit times, and even subdivide catchments into zones with different infiltration and groundwater flow styles, and better understand springshed hydrodynamic behavior.

Acknowledgments The authors thank Miguel Melchor for expert coding and David Mooney and UUCV for their assistance at the Bubble. Funding was provided by the F&M Committee on Grants and GeoScience Founders Society. Data are available publicly (as cited) or from the authors upon request.

Supporting Information Filename Description grl55951-sup-0001-2017GL073790-figs01.docxWord 2007 document , 91.2 KB Supporting Information S1 Please note: The publisher is not responsible for the content or functionality of any supporting information supplied by the authors. Any queries (other than missing content) should be directed to the corresponding author for the article.