Traditional and novel time-series approaches reveal submarine groundwater discharge dynamics under baseline and extreme event conditions

Groundwater is a vital resource for humans and groundwater dependent ecosystems. Coastal aquifers and submarine groundwater discharge (SGD), both influenced by terrestrial and marine forces, are increasingly affected by climate variations and sea-level rise. Despite this, coastal groundwater resources and discharge are frequently poorly constrained, limiting our understanding of aquifer responses to external forces. We apply traditional and novel time-series approaches using an SGD dataset of previously unpublished resolution and duration, to analyze the dependencies between precipitation, groundwater level, and SGD at a model site (Kīholo Bay, Hawaiʻi). Our objectives include (1) determining the relative contribution of SGD drivers over tidal and seasonal periods, (2) establishing temporal relationships and thresholds of processes influencing SGD, and (3) evaluating the impacts of anomalous events, such as tropical storms, on SGD. This analysis reveals, for example, that precipitation is only a dominant influence during wet periods, and otherwise tides and waves dictate the dynamics of SGD. It also provides time lags between intense storm events and higher SGD rates, as well as thresholds for precipitation, wave height and tides affecting SGD. Overall, we demonstrate an approach for modeling a hydrological system while elucidating coastal aquifer and SGD response in unprecedented detail.


Study site
The study was conducted in Wainanaliʻi Pond, located in Kīholo Bay (3.2 km length), on Hawaiʻi Island (19. 27,33 (Fig. 2). Coastal groundwater levels are tidally influenced, even thousands of meters inland from the shore 34 . Within the subsurface numerous structural barriers formed by various types of basaltic lava flows influence, obstruct, and channel groundwater flow, however, these are mostly unspecified 30,35,36 . For example, previous studies have suggested that some Kīholo Aquifer recharge is channeled by a subsurface structural boundary extending from the summit of Hualālai to Puʻu Anahulu 37 while others have confirmed the presence of deep groundwater conduits 31 , yet none of these features and their hydrogeological roles are well described. Nearshore SGD is well-characterized along this coastline, occurring primarily in the form of point discharge (e.g., discrete springs as opposed to diffuse flow) ranging from 200 to 7000 m 3 /day [38][39][40] . Mean annual rainfall (based on the 30-year average) to the watershed ranges from 260 to 1200 mm/year 41 , with the most recharge (288,000 m 3 /day) occurring at higher elevations (500 to 2000 m) along the slopes of the volcanoes 42 . Groundwater pumping (33,

Results
Tidal and seasonal fluctuations. Seasonality in the time series can be broken down into two main oscillations: tidal and seasonal effects. Parameters that were primarily driven by tides were strongly associated with 6-, 12-, and 25-h oscillations and include ocean water level, groundwater level, salinity, wave height, and SGD ( Fig. 3; see Supplementary Figs. S2-S4 online). SGD was greatest at low tide, when the hydraulic gradient (groundwater level minus ocean water level) was at its maximum. Similarly, outside of major precipitation events, salinity was lowest during low tide given the greater influence from SGD. Longer term seasonal effects included air temperature, precipitation, and groundwater level. Wet season (May through October) average air temperature (22 ± 1.2 °C), precipitation (2.1 ± 5.4 mm/hr), and groundwater levels (0.93 ± 0.049 m) were greater than those occurring during the dry season (November through April; 20 ± 0.87 °C, 0.43 ± 1.9 mm/hr, and 0.91 ± 0.029 m, for air temperature, precipitation, and groundwater levels, respectively), although it is worth mentioning that the study period had a significant overlap with an El Niño period. Groundwater level is both tidally and seasonally modulated as it is influenced by both recharge from precipitation and ocean water level. Because of the porous and highly conductive basalt matrix and groundwater elevations near sea level, groundwater levels are affected by ocean water level 34 . While the groundwater monitoring station was 5.6 km inland from the shore and some signal attenuation is anticipated, ocean water and groundwater levels were positively correlated (R 2 = 0.54), with a lag of 5 h (see Supplementary Fig. S5online). The lag is distance dependent and expected to decrease closer to the shoreline. However, neither the distance nor the lag affects the 25-h seasonal and longer-term trends analyses.
Longer term trends. The longer-term trend component was extracted after removing all tidal and seasonal effects (Fig. 4). Here, we focused on interconnectivity between precipitation, groundwater levels, SGD, and salinity. A cursory examination of these parameters plotted with respect to time shows a clear lagged relationship between precipitation and subsequent groundwater levels, offset by a period of two to four weeks. This relationship was confirmed as significant through cross-correlation, with a maximum R 2 value at four weeks of 0.55 (Fig. 5). High groundwater levels were also associated with high SGD three weeks later (Fig. 5).
Relationships between drivers and parameters of interest were further investigated by quantifying the percentage contribution to each parameter using a random forest classifier to calculate feature importance (Fig. 6). Groundwater level drivers included ocean water level, wave height, and precipitation, and then subsequently SGD included all groundwater drivers in addition to groundwater level, and finally coastal salinity analysis included all SGD drivers, in addition to SGD. At the 1-h temporal resolution, tidally driven effects such as ocean water level, wave height, and SGD were the primary influences on groundwater levels. For SGD, tidally influenced variables were again the greatest contributor at the 25-h temporal resolution when El Niño conditions were not present. Conceptual flowchart of processes impacting parameters of interest (indicated in blue boxes -coastal salinity, SGD, and groundwater level). The three parameters analyzed were selected because of their direct and environmentally relevant relationships. Ecosystems are greatly dependent on coastal salinity (parameter 1), which at the study site is driven mostly by SGD (parameter 2) and ocean processes. SGD therefore plays an essential ecological role. Groundwater (parameter 3) is the ultimate source of fresh meteoric groundwater in SGD. Additionally, we must acknowledge that marine drivers are also moderating SGD. Tides are not included as a parameter of interest since tides are predictable and instead are used as an explanatory parameter.

Discussion
Both tidal and seasonal components influence SGD as previously demonstrated in geochemical and hydrological modeling studies 1,3,4,13,14,43 . SGD estimates based off chemical budgets typically represent shorter time scales, whereas hydrological models require adequate knowledge of the aquifer hydrogeology and field observations to create and calibrate the models. Additionally, there are only a few studies that monitor SGD over a sufficient temporal period to discern tidal and seasonal patterns and relationships with its terrestrial and marine driving forces. Here, we propose that SGD observations coupled with precipitation, groundwater level, ocean water level, and wave height data can be used to reveal relationships between SGD driving forces using time series analysis and machine learning methods. Tidal patterns discerned from the time series decomposition (Fig. 4) and lag analysis (Fig. 5) displayed patterns with 6-(M 4 , shallow water overtides of principal lunar tide), 12-(M 2 principal lunar semidiurnal tide), and 25-(K 1 , lunar diurnal tide) hour temporal resolutions. Furthermore, the random forest classification (Fig. 6) demonstrated that short-term tidal effects were still the dominant driver for groundwater levels at the resampled 25-h temporal resolution. This makes sense because groundwater levels are influenced by both tidal (ocean water level Conceptual cross-section of the study area and processes influencing groundwater levels, SGD, and coastal salinity (not to scale). Locations and ages (years before present, BP) of lava flows are approximated on the surface and with depth according to [26][27][28][29] . These historic lava flows are significant as they provide a heterogeneous substrate for groundwater flow that is hard to capture in hydrological models 30 . Precipitation infiltrates into the young, porous basalt and flows in the subsurface before discharging to the coastal ocean as SGD. The assumed location of the subterranean estuary is illustrated by the blue (fresh), purple (brackish), and pink (saline) colors near the discharge zone. These colors are purposefully cut off at a relatively short distance from the shoreline because there is no information on the structure of the aquifer other than the inferred existence of some vertical and horizontal barriers to flow 30 and deeper conduits 31 . The subsurface geology and hydrogeology in the study area is complex and not well-characterized, but it is hypothesized that there are irregularities and/or a confining layer(s) that impede groundwater flow forcing some SGD to discharge 1-2 km offshore (deep SGD; 31 ) in addition to the nearshore, shallow SGD that was the focus of this study 32 . Ocean water level, groundwater level, wave height, tidal pumping, and seawater circulation are all modulated by tidal and seasonal fluctuations. Illustration by Brooks Bays, SOEST Publication Services.  www.nature.com/scientificreports/ and wave height), which are present at 24.8 and 28-day periods, and seasonal (precipitation) oscillations, but as precipitation is episodic in nature, the tidal signature will override the seasonal one. Similarly, tidally influenced parameters had the greatest impact to SGD and salinity during periods that did not coincide with extreme events. Seasonal patterns were associated with fluctuations in precipitation. While precipitation impacts are not instantaneous as demonstrated by the lag analysis (Fig. 5), rainfall, particularly when associated with storm events, led to long-term changes in groundwater levels and SGD (Fig. 4). Over longer time intervals (e.g., comparing 1-h vs. 25-h resolution), precipitation becomes an increasingly important driver to SGD and coastal salinity, as evidenced by the random forest classification (Fig. 6) and studies conducted in other island environments 3 . This corroborates with our lag analyses as precipitation impacts to these components are not instantaneous.
Traditional time series approaches, such as decomposition and lag analysis are useful, but may not be able to fully resolve driver-response relationships between environmental parameters given the influence of overlapping processes. These issues were addressed by conducting additional analyses using a CUSUMs approach to examine driver-response relationships between groundwater levels, SGD, coastal salinity, and precipitation, ocean water levels, and wave height (Fig. 7). The CUSUM plot shows the cumulative sum of the underlying (standardized) time series. A positive CUSUM slope arises from the aggregation of above-average underlying   1,13 . Further evidence was derived from the resulting CUSUM analysis slope, which can be used as a measure of the amplitude of influence in a driver-response relationship. Precipitation rates between 1 and 4 mm/hr for groundwater (1 to 7 mm/hr for SGD) led to the greatest increase in groundwater and SGD rates ( Fig. 7; see Supplementary Table S1 online), where the slope was an order of magnitude greater (72 and 10 × for the relationship between precipitation and groundwater and SGD, respectively) than other precipitation thresholds. The degree of influence declined for precipitation rates between 4 and 34 mm/hr (for groundwater) and 7 and 28 mm/hr (for SGD). Salinity was the lowest at precipitation rates ranging from 1 to 8 mm/hr (slope = − 65x), and then remained below average for precipitation ranging from 8 to 36 mm/hr (slope = 2.2x). The relatively small sample size for precipitation above 40 mm/hr limits our ability to describe driver-response relationships above this threshold. Overall, slope decreased by 96 to 99% from the first precipitation threshold (range of upper bound: 4-8 mm/hr) to the final  Table S1 online). These thresholds indicate the presence of nonlinear response patterns, justifying our multipronged approach. Tidally driven parameters, such as ocean water level and wave height, had a significant impact on groundwater level, SGD, and salinity. Ocean water levels less than − 0.1 m led to low groundwater levels and between -0.1 to + 0.2 m led to high groundwater levels. Ocean water levels greater than + 0.2 m were typically associated with bimonthly spring tides and had a low groundwater levels, suggesting that there is a threshold for tidal modulation of groundwater, perhaps due to constraints on rock permeability 30 . Ocean water levels less than 0 m were associated with high SGD and low salinity, and the opposite was true for ocean water levels greater than 0 m. This is not surprising as SGD is typically greatest at low tide when the hydraulic gradient is the steepest and one would expect salinity to be the lowest when SGD is greatest. Wave heights less than 0.5 m led to low groundwater levels, between 0.5 and 1.3 m led to high groundwater levels, and then greater than 1.3 m led to low groundwater levels before evening out at 1.5 m (Fig. 7). Wave heights between 1.5 and 2.5 m were typically associated with either tropical storms or big wave events. Similarly, while the relationship between wave height and SGD is less clearcut when wave height is less than 0.8 m, wave height between 0.8 and 1.2 m led to high SGD, perhaps reflecting additional influence from saline SGD due to increased wave setup, tidal pumping, and recirculating saline water. Salinity had the opposite relationship to wave height compared to SGD, where wave heights between 0.5 and 0.8 m led to high salinity, while those between 0.8 and 1.5 m led to low salinity. Due to the small sample size for wave heights above 1.5 m, we are not able to interpret above this threshold. The impact of wave height on SGD may sound counterintuitive hydrologically because waves can add water into the coastal aquifer during wave set-up 4 , but if the high waves coincide with spring tides, storms, precipitation events, or higher SGD, then the relationship can be explained by the influence from these processes instead.
Outliers and event detection. Outliers were used for additional analyses and event detection. Univariate outliers (see Supplementary Table S2 online) were examined to look at specific feedback mechanisms but were not used as a basis for outlier detection. Groundwater levels ranged from 0.75 to 1.1 m (average = 0.92 ± 0.053 m with respect to MSL), with outliers occurring for 4 days below the lower bound and 4 days above the upper bound. High groundwater levels occurred during October 2015 and were associated with previous precipitation events sourced from tropical storms in August through September 2015 (see Supplementary Fig. S3 online). Low groundwater levels occurred in May 2014 when ocean water levels were also low. SGD ranged from 0 to 2,000 cm/d (average = 150 ± 240 cm/d), outliers were identified on 3 days below the lower bound and 8 days above the upper bound. Higher SGD rates occurred one to two months after major precipitation events, whereas lower SGD rates were typically associated with lower wave heights, perhaps due to decreased contribution from saline SGD from wave set up. Salinity ranged from 5.8 to 33 during the study period with an average of 18 ± 4.8. Of this, 7 days during the study period were considered outliers below the lower threshold (salinity ranging from 11 to 12) primarily due to high rainfall and 3 were above the upper threshold (salinity ranging from 24 to 25). www.nature.com/scientificreports/ Overall, 25 days were identified as multivariate outliers (Table 1). Of these, 8 days were associated with hurricanes, 12 days were associated with big wave events in February through March 2016, and the remaining were associated with rainfall events. The 2015 Central Pacific Hurricane season was the most active on record, with 16 named storms 6 . On August 30, 2015, a record three Category 4 hurricanes (Hurricanes Kilo, Ignacio, and Jimena) were co-occurring east of the international dateline. While this occurrence was anomalous in the historical record, there is evidence that tropical storms will increase in both frequency and intensity around the Hawaiian Islands 44,45 .
Isolating outliers were informative for extracting relationships that may have an even greater influence in the future. For instance, hurricanes were associated with increased precipitation (in addition to higher seasonal air temperatures and wave height), subsequently leading to increased groundwater levels and SGD (Fig. 3). While ocean water levels typically have an inverse relationship with SGD, this effect is amplified during periods of above average wave height (Fig. 4). This is likely because the hydraulic gradient decreases further, leading to decreased SGD. The effects of ENSO and hurricanes on groundwater levels have been documented in other locations globally, such as the mainland United States and India 46,47 . Similar to our results, a model-based study reported SGD, groundwater level, and recharge anomalies that were linked to ENSO for a barrier island in North Carolina, USA 48 .
Climate change in Hawaiʻi is anticipated to lead to decreased monthly average precipitation (but an increase in extreme rainfall events), increasing sea levels, increasing water and air temperature, disruptions to trade wind patterns, and an increase in the hurricane frequency locally 44 . The data collected during this study provides vital information about how precipitation can impact groundwater levels, SGD, and coastal salinity. For instance, an increase in the frequency and intensity of catastrophic rainfall events (and coincident decline in average Figure 7. CUSUM driver-response plots (columns = drivers, rows = response). Drivers shown include precipitation rates exceeding 1 mm/hr (mm/hr; column 1), ocean water level (OWL; m with respect to MSL; column 2), and wave height (Wave Ht; m; column 3). These are compared to the following responses: groundwater levels (m; row 1), SGD (cm/d; row 2), and salinity (row 3). The underlying data showing the relationship between x and y are shown in the green dots on the primary y-axis. The driver-response relationship is illustrated with the black line (secondary y-axis), comparing the CUSUM groundwater level, SGD, and salinity vs. precipitation, OWL, and wave height. www.nature.com/scientificreports/ precipitation) could lead to seasonal changes in the availability of water resources or the magnitude of SGD. All recorded tropical storms that passed near the main Hawaiian Islands during El Niño were associated with outliers in this study, particularly due to changes in precipitation and groundwater levels. Multiple studies have demonstrated that decreased precipitation and increased stress on water resources leads to decreasing fresh SGD, higher sea levels, and salinization of coastal aquifers 49,50 . Major rainfall events also lead to decreasing coastal salinity from storm runoff, which can have negative ecological impacts as biota tend to be sensitive to certain salinities and sediment load. Both in Hawaiʻi as well as globally, a decrease in fresh SGD and increasing sea levels could have major negative impacts to coastal ecosystem health. Groundwater dependent ecosystems, which rely on lower salinity SGD with optimal nutrient levels, are aggravated by changing precipitation patterns and sea-level rise 51 . Humans have relied on these ecosystems for years-for instance, traditional Hawaiian aquaculture depended on SGD to feed coastal fishponds 52 , such as the one in Kīholo Bay. Beyond looking at this from a human-impact perspective, anthropogenic influence on SGD has also led to coastal eutrophication and harmful algal blooms in for instance, Bermuda 53 , and has been associated with an increase in coral reef disease severity in Guam 54 . These groundwater dependent ecosystems are vulnerable and are posed to face additional challenges in the upcoming decades.
The methods and analyses used here are applicable to other coastal aquifers globally but are particularly useful for areas that lack well-characterized hydrogeology and thus are difficult to describe with a hydrological model. Here, we were able to constrain the temporal relationships between variables, elucidate driver-response relationships, determine which predictors have the greatest influence on a variable, and isolate events to learn more about the relationship between atmospheric, terrestrial, and oceanic processes under conditions that will likely be more common in the future. In this way, our methods provide similar results to a hydrological model without the requirement of a highly monitored subsurface. Given the conditions of the field study, we were unable to separate the fresh and saline SGD fractions due to variations in the radon to salinity ratio and the embayment residence time over longer time periods. While the radon to salinity ratio holds consistent for shorter time periods (weeks to months), seasonal factors such as rainfall or large swell events make it challenging to derive the fresh SGD fraction without very high uncertainties 55 . In another example, urban coastal aquifers frequently experience a "coastal groundwater squeeze" sourced from declining water resources due to increased salinization, overextraction, land subsidence, and pollution 2 . Additional applications could include providing insight to feedback relationships between coastal salinity and precipitation and/or sea-level rise in terms of fish abundance, supply, and diversity. In all, the methodology applied in this study can provide valuable predictive feedback for understanding perturbations to the water cycle in a wide variety of environments.

Conclusion
This study applies both traditional (e.g., decomposition, cross-correlation) and novel (e.g., random forest classification, CUSUMs analysis) time-series approaches to understand coastal aquifer responses to atmospheric, terrestrial, and oceanic forces at a model study site with poorly characterized hydrogeology. Using these techniques, we were able to establish that there is a two to four week temporal lag between precipitation and groundwater levels. Increases in SGD occurred three weeks following higher groundwater levels. We were also able to demonstrate that tropical storms led to increased SGD due to increased precipitation. The frequency and magnitude of these tropical storms also coincided with a strong El Niño, and the influence of precipitation on SGD increased by over 9% during this period. Yet, outside of anomalous events, tidally driven parameters had the greatest instantaneous influence on groundwater levels, SGD, and coastal salinity.
This research demonstrates the utility of combining field-based hydrological measurements with data-driven approaches, particularly for areas where the underlying aquifer characteristics are difficult to portray in a hydrological model. The methodologies applied in this research can be used to provide further insight into water

Methods
Data collection and SGD calculation. Radon data were collected using an autonomous gamma spectrometer called the SGD Sniffer at 0.3 m depth below the surface of the water 32 , which was deployed from March 2014 through July 2016 in Wainanaliʻi Pond. Conductivity, water depth, and water temperature were concurrently logged with a CTD probe (Schlumberger Inc. CTD diver). All data were collected with a 1-h time resolution. Briefly, SGD was calculated using a transient radon mass balance (Eq. (1); see Supplementary Fig. S6 online 56 ), where F SGD is the SGD advection rate (m/day), A is the radon activity (Bq/m 3 ) measured by the SGD Sniffer, A gw , A ocn , and A atm refer to radon activities in groundwater, offshore ocean, and the atmosphere, respectively, A Ra is the dissolved 226 Ra activity, z pl is the plume thickness (m), T refers to the water residence time (days), and F atm and F mix account for radon losses due to atmospheric evasion and lateral mixing (Bq/m 2 /day). Uncertainties were propagated throughout the mass balance calculations.
F atm was calculated using Eq. (2), where k is the gas transfer coefficient (m/day), a is the radon gas solubility coefficient, and A atm is the radon activity in the atmosphere. Local wind speed (NOAA Tides and Currents, m/s) was used to derive the radon gas transfer and solubility coefficients.  Supplementary Table S3 online). To overcome these periods of missing data, we primarily focused on using time series analysis methods that can accommodate missing data.
The study period partially overlapped with the 2014-2016 El Niño, with an ONI ranging from 0.6 to 2.6 6 . Additionally, the study period coincided with the 2015 Pacific hurricane season, the second most active on record 6 and the most active within the central Pacific 6 . In winter of 2016, big wave events also overlapped with the monitoring period.
Data pre-processing. All analyses were conducted in Python and R. First, data were power transformed by adding 1 and taking the square root so that its distribution approximates the normal distribution more closely. This was done for all analyses except the cumulative sums (CUSUMs) analysis (described below). After this, data were standardized using Eq. (3), where x i refers to the i-th value in the time series, m is the mean, σ is the standard deviation, and Z i represents the standardized value.
Because hydrological processes and interactions may be more apparent on different time scales or are not instantaneous events, we also conducted a second set of analyses using data that were resampled to a 25-h temporal resolution with short-term tidal fluctuations removed. The tidal signal was removed from the 25-h resampled data by calculating a running average over a six-hour interval and subsequently a 25-h interval. The data were then resampled into 25-h time intervals using the median value from each period.
Time-series analyses. Following initial preprocessing of data, we conducted an additive time series decomposition using moving averages for each variable using statsmodels 60 to derive the trend, seasonal, and residual components. From here, we extracted the de-seasoned trend component for further analyses to look at longer-term relationships. The Augmented Dickey-Fuller Test was used to confirm stationarity of the decomposed components at the 5% level of significance.
Outliers were detected using two different methods-univariate outliers were identified using the 1.5 interquartile range (IQR) rule and multivariate outliers were detected by calculating the Mahalanobis Distance. The univariate outliers were not removed from the data, but instead were used to observe individual predictor variance and trends. The multivariate outliers were removed from the data and set aside for additional event classification and analysis. In most cases, however, there was substantial overlap between outlier date and time identified using either method. To detect univariate outliers, the 1.5 IQR rule was used, where a data point is considered an outlier if it is less than Q1 -1.5 * IQR or greater than Q3 + 1.5 * IQR, where Q1 and Q3 refer to the first and third quantile, respectively.
Multivariate outliers were detected by calculating the Mahalanobis Distance with the de-seasoned dataset using Eq. (4) (MD; 61 , where the multivariate vector is represented by x , with mean µ , and S is the covariance matrix. For each time point, p-values were calculated for the resulting MD by calculating the χ 2 value for each MD. Outliers were identified when the resulting p-value was less than 0.001. www.nature.com/scientificreports/ Temporal lags between precipitation, groundwater, and SGD were assessed for the 25-h resampled data using the cross-correlation function (CCF) in R. Cross-correlation plots can indicate either lagged or instantaneous relationships between two variables. The estimated lag time was determined as the lag with the greatest R 2 .
To establish the relative influence of each driver, we conducted a random forest classification to extract feature importance for groundwater levels, SGD, and coastal salinity. This analysis was conducted using the scikit-learn RandomForestRegressor 62 in Python for both 1-h resolution and 25-h resampled data. Feature importance was examined during El Niño conditions by constructing separate random forest classifiers for data occurring when the ONI was greater than 0.5 (indicating El Niño conditions) and less than 0.5 (indicating El Niño conditions were not present).
CUSUMs analysis was used to evaluate driver-response relationships. This method can help resolve variable relationships in complex environmental time series data with missing intervals 63 . To quantify driver-response relationships, data were first standardized using Eq. (2). After data were standardized, driver (precipitation, ocean water level, and wave height) and response (groundwater levels, SGD, and salinity) relationships were quantified. For each driver-response relationship (e.g., precipitation impact on groundwater levels), driver data were first organized in ascending order. Then the CUSUM for the ith observation for the reordered response variable, s i , was calculated using Eq. (5).
The standardized data, z i , have zero mean and unit variance (m = 0, σ = 1). When the z i value is positive, the underlying value, x i , is greater than the mean. Similarly, when z i is negative, then x i is less than the mean. The slope of the CUSUM trend is an indication of the sign and magnitude of the x i values: a steep positive slope is the result of large positive deviations of x i from its mean, and a steep negative slope is the result of large negative deviations of x i from its mean. CUSUM analysis can be used to examine relationships between variables 63 .