Recent groundwater and lake-stage trends in Cape Cod National Seashore: relationships with sea level rise, precipitation, and air temperature

Hydrology influences numerous physical, chemical, and biological processes in aquatic systems, (Leira & Cantonati 2008). Seasonal fluctuations in water depths occur in virtually all natural lakes and wetlands, and are a part of the normal functioning of these systems. Long-term, multi-year or -decadal trends in water levels, however, may be the result of anthropogenic influences, including climate change (Bates 2009, and references therein). Across much of the northeastern United States, the amount of precipitation and frequency of extreme events has increased (Groisman et al. 2001; Huntington et al. 2009). Moreover, rainfall is predicted to increase in this region by between 7 and 14% annually by 2099, particularly during the winter (Hayhoe et al. 2007, 2008; Dukes et al. 2009). In addition, responses of water levels to precipitation vary with numerous factors, including seasonal variability in rainfall, lake morphometry, runoff, inflow/outflow (Vassiljev 1998), and evaporation/evapotranspiration (especially drought), both of which are increasing in many regions with atmospheric warming (Adrian et al. 2009).

Coastal aquifers are also influenced by sea level rise (SLR). Fresh groundwater (GW), which sits atop the denser seawater, is pushed upward as sea level rises – a response that has been observed in many regions around the world (Noest 1991; Turner et al. 1997; Singh 2002; Mohal et al. 2006, Shamsudduha et al. 2009; Bjerklie et al. 2012; Rotzoll & Fletcher 2013; Dong et al. 2015; Hoover et al. 2017). Even short-term events such as wave run-up (a series of successive high waves) and Spring high tides can cause groundwater to rise and stay elevated for a period of time (Lanyon et al. 1982; Hegge & Masselink 1991; Campbell & Bate 1998; Rotzoll et al. 2008). According to the Intergovernmental Panel on Climate Change (IPCC), global sea level is predicted to rise between 52 and 98 cm by the year 2100 which translates to a possible doubling or tripling of the long-term rate (1965–present) of ∼3.65 in this region to ∼6–11 mm/yr (IPCC 2013). Of course, SLR effects on groundwater are spatially variable. For example, Song & Zemansky (2013) reported that inland groundwater levels are affected less than those near the coast. Groundwater modeling of mid-Atlantic coastal systems showed that SLR will elevate GW levels by about 50% of SLR (Wang et al. 2017).

Another important factor influencing groundwater and pond water levels is evaporation/evapotranspiration, processes which are increasing with a warming climate (Szilagyi et al. 2001; Abtew & Melesse 2013; Smith et al. 2016). The inverse relationship between water levels and temperature tends to be stronger during periods of hotter weather (Chen et al. 2004). Also, temperature can influence groundwater levels much more than precipitation where the depth to aquifer is shallow (Chen et al. 2004). Thus, higher air temperatures may increasingly offset groundwater recharge and effectively lower groundwater levels if precipitation patterns remain relatively constant.

Cape Cod is a coastal peninsula in southeastern Massachusetts (USA) with hundreds of freshwater lakes. Formed by sediment runoff from the Laurentide glacier about 23,000 years ago, most of the landscape is classified as ‘outwash plain’ and consists primarily of sand and gravel. On an area basis, the six unconfined groundwater aquifers and surface lakes are the most important hydrologic features on Cape Cod (Figure 1). The latter, originally formed from depressions in the landscape created by large blocks of melting ice left behind by the retreating glacier, have roughly circular shapes and are known locally as ‘kettle ponds’. Many hundreds of these ponds occur throughout the peninsula and they are highly valued as aesthetic, cultural, and natural resources. Aside from the recreational opportunities they provide, kettle ponds are a critical habitat for various species of flora and fauna. Several aquatic plants designated by the State of Massachusetts as Endangered, Threatened, Special Concern, and Watch Listed are found along their shorelines (LeBlond 1989). The ponds also support a number of state-listed odonates and gastropods as well as more common species of reptiles, amphibians, fish, birds, and mammals.

Given that Cape Cod lacks any significant riverine inputs, rainfall is the sole source of freshwater on Cape Cod and the dominant factor influencing Cape Cod National Seashore (CCNS) kettle pond and groundwater levels (Craine & Orians 2004; Masterson & Portnoy 2005; Massey et al. 2006), However, increases of as little as 20 cm in sea-level can have large effects on groundwater according to Masterson et al. (2014). Based on observations from a United States Geological Survey (USGS) well within CCNS, groundwater rose ∼2.1 mm/year between 1950 and 2000 (McCobb & Weiskel 2003). Walter et al. (2016) modeled groundwater responses to SLR within the Sagamore and Monomoy lenses on upper Cape Cod, and suggested that the increases in the groundwater table resulting from 6 feet (182.88 cm) of SLR would range from almost no response to nearly 6 feet, depending on proximity to the coast and to various surface water features.

On outer (also known as ‘lower’) Cape Cod, both groundwater and pond water levels have been measured within CCNS on a monthly basis from 2000–2017 (Figure 1). In this study, we analyzed temporal trends in eleven groundwater wells and nine kettle ponds within CCNS during this period. We further explored relationships between the hydrology of these systems with local precipitation, air temperature, and SLR during this period. The results are important for (1) understanding how CCNS's aquatic resources are changing in a physical sense, (2) predicting how such changes may influence the ecological functioning of these resources, and (3) developing appropriate management strategies, if feasible.

Groundwater elevations were measured using an electric tape that produces an audio signal when the probe contacts the water surface upon being lowered into the well. Depth to water was recorded and then water levels were converted to a known elevation (NAVD88) based on nearby benchmarks. The wells were monitored monthly and the data from all ponds collected within the same 24 hour period. The wells were not sampled within 24 hours of a rain event. The most complete period of record spanned the months of April through November as there were a number of missing data during the months of December, January, February, and March due to ice cover. One USGS well (TSW-89) situated within the park was included in the analysis since there were data spanning the same period of time within each year (April–November) between 2000 and 2017. The data for this well are presented as GW elevations relative to the lowest level recorded during that period of time since the raw data are in feet below ground surface. Unfortunately, a number of other USGS wells within the park boundary could not be analyzed due to major differences in the frequency and timing of sampling, as well as a certain amount of missing data. This prevented direct comparisons of most USGS wells with CCNS wells.

For CCNS pond stage monitoring, siphon gauges were used from 2000 through 2005 (McCobb et al. 1999). A siphon gauge consists of a standpipe on the shore of the pond that connects with a flexible hose to a point in the waterbody. The water level in an onshore well and the surface water level of the pond equilibrate once a connection is made and a reading using an electric-tape method is taken. In 2005, a new method was implemented based on differential leveling. An auto-level was set up between the water body and the existing siphon gauge (which served as a benchmark elevation). Differential leveling was then used to correct pond stages to known elevations.

While analyzing the data, it became apparent that plots of continuous, high-resolution groundwater elevations from well WNW-124 and Great-T pond stages in relation to precipitation during the Fall–Winter of 2009–2010 could shed some light on the sensitivity of the former to the latter. During that time, water levels were recorded using HOBO™ pressure-loggers programmed to collect data every 15 minutes, which were then barometrically corrected based on an additional unit left in the air.

NOAA tidal data were acquired from a tide gauge in Nantucket Harbor (the nearest gauge to the CCNS kettle ponds) for 2000–2017 (https://tidesandcurrents.noaa.gov). This dataset consisted of monthly values in m NAVD88 (the North American Vertical Datum of 1988) for the mean highest high water (MHHW), mean high water (MHW), mean sea level (MSL), mean low water (MLW), and mean lowest low water (MLLW). Precipitation and air temperature (daily maximums, averages and minimums) datasets for the 2000–2017 period were obtained from a weather station located in Chatham, approximately 25 km away from the cluster of CCNS kettle ponds (Station KCQX, https://www.ncdc.noaa.gov). These data were analyzed for temporal trends over various time periods described below. The precipitation data includes snowfall, which was converted to millimetres of water for analysis. Precipitation data represent rainfall totals during the specified period, whereas mean values for tide heights and temperatures were calculated and analyzed as described below.

Precipitation totals were calculated for various time intervals to examine seasonal patterns. These included all data (from May of the previous year through April of the current year; or ‘annual’), November–April (main recharge period), November–August (recharge + drawdown period), June–August (summer; period of maximum biological activity), and August alone (period of lowest water levels). The rationale for these divisions was that rainfall during certain times of the year could influence water levels more than others. November–April is generally when most recharge from rainfall occurs and evaporation and evapotranspiration are negligible (Masterson & Portnoy 2005). Frequencies of precipitation events >0, 10, 20, and 40 mm were calculated for each year, as was the duration of rainfall events in days.

Finally, the total changes in GW elevation change during 2000–2017 were plotted against distances of the wells from both the coastline and the center of their corresponding GW lenses. The latter two variables were calculated using the measuring tools in ARGIS ver. 10.

There were a number of years when GW elevations and pond stages were not documented in certain months, largely due to freezing temperatures/ice cover in December through March. On dates where a number of wells or ponds were sampled, but not others, missing values were replaced with predicted values based on wells that are highly correlated with each other (R² > 0.80). Accordingly, missing values for one well could be estimated where there were data for the other well in a specific pairing, and vice versa, according to the regression equations developed for specific well or pond pairings. For dates when no wells were sampled in between two months, the averages of the previous and successive months were calculated (i.e., mid-points). Overall, 10.3% and 4.5% of the groundwater and pond stage datasets, respectively, were comprised of these predicted values (a list of missing data that was replaced in this way is presented in Table 1).

For analysis, mean values for each variable were calculated by year (and segments of years such as April–September) to reduce seasonality effects. The population means of all variables were all normal and heteroscedastic as revealed by distribution and homogeneity of variances tests (JMP ver. 10.0.2). This allowed for linear regression (α = 0.05) against time and the other variables. Mixed-model stepwise regressions (with p values to enter/leave in the equations set to 0.25) were used to ascertain how precipitation, sea level, and temperature contributed collectively to GW and pond stage variability and trends. The three variables (components of tide heights, precipitation, and temperature) that demonstrated the best individual fits with GW or pond stage were used in the model. Loess curves of water level, precipitation, and tidal data were calculated and plotted together to illustrate relationships among these variables (JMP ver. 10.02).

While R² values were low, GW elevations in 4 out of 11 wells (TSW-106, TSW-134, TSW-89, and WNW-122) exhibited statistically significant increases between 2000 and 2017 and the rest, while statistically insignificant, showed increasing trends (Figure 2). Based on the regression equations, the average magnitude of GW elevation change (mean of all wells) was 8.2 mm/year, ranging between 2 mm/year (TSW-258) and 11.3 mm/year (TSW-89).

No ponds exhibited statistically significant trends during the time period analyzed, even though six of nine ponds exhibited apparent increases (Figure 3). In contrast, Duck, Gull, and Herring ponds declined slightly (although not significantly so) between 2000 and 2017.

Average precipitation by month throughout 2000–2017 was highest in March and November, and lowest in June and July (Figure 4). None of the precipitation totals for the various time periods analyzed exhibited significant trends between 2000 and 2017. Mean precipitation amounts (excluding zero values), however, did increase significantly, while the average duration of precipitation declined (Figure 5). The frequency of precipitation of any amount also declined during this period but the frequency of more extreme rainfall events (>10 to >40 mm) demonstrated no significant trends (Figure 6).

Since 2000, MSL has risen by approximately 0.13 m (= 6.7 mm/year), although the long-term rate since 1965 is 3.65 mm/year (https://tidesandcurrents.noaa.gov/sltrends) (Figure 7). The recent rise above the long-term rate is presumably because the rate of SLR has increased in recent decades but may also be part of larger decadal-scale fluctuations in sea level as well (Church & White 2006; Nerem et al. 2006). Regardless, the magnitude of increase in tide heights between 2000 and 2017 was only slightly variable among the different tidal parameters, ranging between 0.09 m (5.0 mm/year) for MLW and 0.14 m (7.8 mm/year) for MHW (Figure 7).

When hydrographs of well WNW-124 and Great-T Pond were plotted with precipitation, the latter responded quickly and noticeably to rainfall events >20 mm, but not to more typical rainfalls that were below this threshold (Figure 8). In contrast groundwater level exhibited almost no discernable response to discrete rainfall events. Instead, groundwater appeared to be reacting to cumulative rainfall over the course of months. This is demonstrated by the upward trend of groundwater through multiple rain events throughout several months (Figure 8).

Groundwater elevations were best correlated (positively) with November–April MHHW, except for TSW-257, which only exhibited a significant (positive) relationship with November–April total precipitation (PPT). Correlations between GW levels and November–April average daily temperatures were generally very low and insignificant, with the exception of TSW-106 and TSW-89, both of which curiously exhibited a statistically significant positive relationship. For every well, however, combinations of variables, particularly MHHW and PPT, in mixed-model multiple regression analyses produced higher R² values (Table 2). Figure 9 illustrates similarities and differences in temporal patterns of change, plotted as Loess-smoothed curves, among TSW-106, PPT and MHHW (temperature was not included due to this variable generally being excluded by multiple regression models). Of particular note is the way in which GW fluctuations generally track changes in MHHW through time.

With the exception of Great-W, Gull, and Herring ponds, water levels were significantly correlated with November–April PPT (Table 2). Only Dyer and Great-W demonstrated significant relationships with November–April MHHW (positive) and no ponds exhibited relationships with November–April average temperature. In general, combination of PPT and MHHW yielded much higher R² values compared with individual variables (Table 2).

When the total change in GW elevations over the 2000–2017 time period calculated from regression equations (statistically significant or not) was plotted against distance from the coastline, there were no significant relationships (data not shown). However, when these data were regressed against distance from the center of the GW lens within which the wells reside, a significant positive relationship emerged, suggesting that wells furthest away from the center of the lens (top of the aquifer; highest elevation of GW within the substrate) exhibited a larger rise than those nearer the center (Figure 10). Two wells (WNW-105 and WNW-122) were excluded from this analysis since they are both sandwiched between the Pamet and Chequesset lenses, and are not exclusively part of either system (Figure 1).

Groundwater within CCNS has risen significantly in several wells within CCNS, despite no detectable increases in precipitation between 2000 and 2017. Wells demonstrating significant upward trends (TSW-105, TSW-134, TSW-89, and WNW-122) were the first, second, third and fifth closest to the coastline, respectively, out of the eleven wells. Moreover, there was a significant, positive relationship between the magnitude of GW rise and distance from the center of the corresponding GW lenses, where subsurface water elevations are highest (Figure 10). This generally agrees with the modeling of Walter et al. (2016), which suggests that GW closer to the coastline (i.e., further away from the top of the aquifer) will be more responsive to sea level change rise than sites further inland and/or nearer kettle ponds and surface water outlets (rivers, streams).

Although R² values for temporal trends in GW elevations were very low, this can be expected given the influence of so many other factors. Surface/subsurface geology, topography, land use, groundwater extraction, transmissibility, water storage capacity, specific yield, position on the GW lens, hydraulic gradient, depth to the groundwater table, and especially seasonal variability can all contribute to substantial variability around the general trend (Gillham 1984; McCobb & Weiskel 2003). However, low R² values can still indicate a genuine relationship between variables (Rubinfeld 2000). Moreover, every GW well demonstrated a positive change, irrespective of statistical significance. Based the consistency of the increasing responses among wells, it appears that GW elevations in CCNS may be at least partially responding to sea level change. Between 2000 and 2017, MHHW, MHW, and MSL rose 0.11 m, 0.14 m, 0.11 m, 0.09 m, and 0.10 m respectively (= 6.1 mm/year for MSL), although the longer-term trend is much lower (3.57 mm/year between 1965 to present; https://tidesandcurrents.noaa.gov/sltrends). For reference, the average rise in CCNS GW wells was 6.3 mm/year.

That GW trends were not offset to a larger extent by a warming climate during this time, which enhances evaporation/evapotranspiration, is somewhat surprising. Although summertime air temperatures have risen ∼2 °C (Smith et al. 2016; Smith 2017) since 1983 and PPT has remained relatively stable over the time period analyzed, WNW-122, TSW-105, and TSW-134 all exhibited significant increases in GW elevations during the month of August (data not shown). In addition, the inter-annual variability of the GW data was amplified by the very low values in 2002, when the most severe drought in this time period occurred (https://www.drought.gov), and the following year in which precipitation and GW were extremely high. In addition, diminished wind speeds (due to growth of the surrounding forests and/or changing weather patterns) results in lower rates of evaporation during the warmer months (Pilson 2008; Smith et al. 2018).

There are circumstances under which groundwater can rise at the same rate as sea level (Manda et al. 2015; Hoover et al. 2017). Rozell & Wong (2010) modelled a 0.61 m sea-level rise with a simulated 2% reduction in precipitation for a sandy island in Long Island Sound and reported a subsequent increase in GW elevation of 0.59 m. This was partially explained by deeper clay layers that limit aquifer depth, which may be relevant to CCNS as well given that Portnoy et al. (2001) described many ponds on the outer Cape as being perched to some degree due to impermeable clay layers below them. Oldale & Barlow (1986) and Zeigler et al. (1964) also reported that portions of the Wellfleet plain deposits contain beds of clay. Walter et al. (2016) modeled the responses of upper Cape Cod aquifers to SLR scenarios of 0.6–1.8 m and found that the median rise in GW across the Sagamore lens was roughly 25–30% of SLR. Within the much smaller Monomoy lens, however, median values for GW rise were estimated to be 45–50% of SLR, while areas adjacent, or in close proximity, to the coast would experience nearly the same rate of rise as sea level. Higher sensitivities of small aquifers to SLR have been reported previously for other coastal regions (Chui & Terry 2013). Given that the Chequesset and Pamet lenses on the outer Cape (within CCNS) are nearly five times smaller in area than the Sagamore lens of the upper peninsulas, we might therefore expect the groundwater elevations in this area to be more responsive to sea level in terms of magnitude of change. Unfortunately, because of the timing and frequency of data collection for the USGS wells in this region of the Cape they could not be directly compared to CCNS well data.

There was also likely a much larger effect of rivers and streams that drain water out of the kettle ponds across upper and middle Cape Cod. These features, which significantly dampen GW responses to SLR (Masterson & Garabedian 2007; Walter & Masterson 2011), are largely absent on the outer Cape. The hydrologically restricted Herring and Pamet Rivers are the only surface water outlets for the Chequesset and Pamet lenses, respectively. The land surface over of the Sagamore and Monomoy lenses on upper and middle Cape Cod also is more highly developed with extensive impervious surfaces that enhance direct overland runoff of precipitation rather than percolation through the substrate. In addition, the clear-cutting of forests and woodlands for urban development tends to result in more groundwater evaporation (Taylor & Stefan 2009). Finally, there is far more groundwater extraction on the upper vs. lower Cape for public use.

Spatial variability in rainfall can be considerable (Chaubey et al. 1999; Tetzlaff & Uhlenbrook 2005; Tokay et al. 2014), particularly in coastal vs. inland areas (Daly 2006; Kafle & Bruins 2009). In CCNS, a decrease in the frequency of rainfall was observed although average precipitation amounts have actually increased. The frequency of large rainfall events (i.e., 10–40 mm) either declined or exhibited no change within the time periods analyzed. This contrasts with patterns observed throughout many parts of the northeastern US, where extreme rainfall is becoming more common (Douglas & Fairbank 2010). This is important because larger downpours contribute most to groundwater recharge in many instances where there is minimal runoff (Nimmo et al. 2005; Owor et al. 2009; Taylor et al. 2013; Zhang et al. 2016). For CCNS, the continuous water level data indicated that only rainfall events of >10 mm or more perceptibly impact pond water levels in the short term, and there are no immediate changes in GW elevations after such events. Antecedent moisture conditions and the thickness and nature of the unsaturated zone, including the capillary fringe, are also major factors affecting recharge (Gillham 1984).

According to Walter et al. (2016), certain hydrologic features such as kettle ponds and rivers/streams serve as ‘relief valves’ for rising GW tables, which may explain their lack of significant rising trends. Herring and Gull pond stages, which exhibited significant reductions in GW elevations between 2000 and 2017, appear to reflect this process as they both drain into the Herring River (Figure 1). Rising groundwater around these ponds would theoretically cause more water to enter the river and be discharged to the coast (Masterson & Garabedian 2007). The nearest groundwater wells (WNW-105, WNW-122) to these ponds exhibited very weak correlations with pond stages – unlike the pairings of other wells in close proximity to individual ponds that are ‘landlocked’. They are also situated at the convergence of the Pamet and Chequesset lenses, where GW elevations are lowest and where groundwater extraction occurs at a town (Wellfield) well field (Figure 1).

Fordyce (2014) suggested that SLR best explained the long-term rise in groundwater, while rainfall events accounted more for short-term fluctuations in a New Zealand coastal aquifer. The results of this study similarly suggest that decadal-scale increases in certain groundwater wells, may at least partly be driven by SLR, while the seasonal and inter-annual variability of precipitation contributes to short-term fluctuations, or ‘noise’ around the general trends (Figure 9).

Regardless of the mechanism(s) by which ground- and pond water levels have changed between 2000 and 2017, altered hydrology may be having important shorter-term ecological consequences. For example, if a pond shore littoral zone has an elevation gradient of 5%, a narrowing of 2.3 m of this zone (horizontal distance) would be expected from a water level rise of 20 cm. Prolonged high water levels during a single season can reduce aquatic macrophyte biomass and species richness in the littoral zone (Farney & Bookhout 1982; Sjöberg & Danell 1983; Ashton & Bissell 1987; Nohara 1991; Casanova & Brock 2000; Geest et al. 2005; Trei & Pedusaar 2006). In addition, diminishing drawdown in summer can adversely impact species that require exposure for seed germination and early development (Keddy & Reznicek 1986; Wagner & Falter 2002; Leira & Cantonati 2008). On a landscape scale, rising groundwater elevations (if these trends continue) will eventually result in the development of new wetlands in topographic depressions that are currently dry. Changes in hydrology of extant seasonal wetlands can also be expected – some of which will eventually convert to permanent ponds. Moreover, elevated GW generally equates with increased discharges through streams and rivers (Nuttle & Portnoy 1992), which could then impact coastal water quality due to elevated nutrient inputs from fertilizers and septic systems (Bowen & Valiela 2001; Hauxwell et al. 2003).

Conversely, the hydraulic head differential between the groundwater table and sea level diminishes with increases in the latter (Holliday et al. 2007; Anderson 2015; Gonneea et al. 2013). Thus, reduced groundwater flow through the substrate is expected with SLR, which can subsequently affect coastal marshes, seagrasses, and estuaries by altering nutrient inputs and salinities (Johannes 1980; Nuttle & Portnoy 1992; Rutkowski et al. 1999; Bowen et al. 2007; Masciopinto & Liso 2016). This may occur to a greater extent on the outer Cape (e.g., within CCNS) where the smaller Chequesset and Pamet aquifers seem to be more sensitive to sea-level perturbation. Of course, all these relationships could change if patterns and amounts of precipitation change with evolving climate conditions.

In conclusion, while the hydrologic trends observed in this study are fairly subtle and the result of complex processes within the substrate, groundwater has generally risen over the last 17 years (significantly so in locations closer to the coastline). SLR during this time may be partly responsible for the longer-term trends, while precipitation contributes more to variation about these trends. On the other hand, pond stages have remained relatively constant and appear to be driven more by direct precipitation. Temperature increases apparently played very little role in regulating water levels during 2000 to 2017. However, that could change with continued climate warming. Higher temperatures have already been shown to influence kettle pond thermal structure (Smith et al. 2016, 2018). For CCNS management, such processes are beyond the ability of the park to manage. However, the results serve as a basis with which to interpret any changes in biological variables (e.g., plant species composition and/or densities) that have occurred during this time, and how these ecosystems may respond to further accelerated SLR in the near future.

This work was supported by the Department of the Interior, National Park Service, Cape Cod National Seashore.