Spatially-Resolved Integrated Precipitation-Surface-Groundwater Water Isotope Mapping From Crowd Sourcing: Toward Understanding Water Cycling Across a Post-glacial Landscape

Isotopic analyses of δ18O and δ2H of water in the context of the hydrologic cycle have allowed hydrologists to better understand the portioning of water between the different water domains. Isoscapes on a large spatial scale have been created to show isotopic variation in waters as a function of elevation, temperature, distance to coast, and water vapor source. We present the spatial and temporal isotopic results of precipitation, surface water, and groundwater of an ongoing study across Massachusetts, USA in order to establish an isotopic baseline for the region. This represents one of the most comprehensive and detailed isotopic studies of water across a 10,000 sq mi area that has exhaustively sampled important components of the terrestrial hydrologic cycle (precipitation, groundwater, and surface waters). We leverage the support of volunteers and citizen scientists to crowd source samples for isotopic analysis. The database consists of water samples from 14 precipitation sites, 409 ground water sites and 516 surface water sites across the state of Massachusetts, USA. The results indicate that groundwater isotopic composition ranges from δ18O −11 to −4‰ surface water ranges from δ18O −13 to −3.84‰ and precipitation ranges from δ18O −17.88 to −2.89‰. On a first order, the small bias of mean groundwater (−8.7‰) and surface water (−8.0‰) isotopes compared to precipitation δ18O (−7.6‰) supports that groundwater recharge and surface water storage effects through the hydrologic year impact the isotopic composition of surface and groundwater. While differences are distinct, they are larger than previously reported values, but still suggest more importance of summer precipitation than previously acknowledged. On average seasonal amplitudes of precipitation (2.7‰), surface water (1.13‰), and groundwater (~0‰) of the region demonstrate young water fractions of surface water to be 40% with groundwater ~0%. Results demonstrate that mean δ18O in precipitation, surface water and groundwaters are more enriched in heavy isotopes in areas near the coast, than the interior and western portion of Massachusetts. The hope is for this dataset to become an important tool for water management and water resource assessment across the region.

Through the creation of isoscapes, defined here as a spatial and temporal distribution of water isotopes for a particular geographic region, local processes such as the interaction between surface and groundwater can provide spatial and temporal information that can be beneficial to water resources and management. Isoscapes are the end result of spatial and temporal distribution of isotopes, they are useful in describing environmental conditions across space and time (Bowen, 2010). They allow us to determine the interconnectivity in various systems, such ecological, hydrological, biogeochemical, or meteorological systems. Kendall and Coplen (2001) created an isoscape for δ 18 O and δ 2 H in river waters across the United States using precipitation isotope data collected from selected US Geological Survey water-quality monitoring sites. Through extensive analyses of δ 18 O and δ 2 H in both precipitation and surface waters this study revealed distinct spatial and temporal differences in δ 18 O and δ 2 H of both river and meteoric waters. State Meteoric Water Lines were created based on geographic regions across the US. A regional pattern within the State Meteoric Water Lines was noted suggesting the large geographic areas are controlled by the humidity of the local air mass which conveys an evaporative enrichment and thus in the stream samples within the area. Spatially the δ 18 O and δ 2 H in river waters fall along local water lines but also reveal a distinct correlation with topography and appear to be primarily reflecting the isotopic signal of precipitation. Jasechko et al. (2014) used 20 years of spatial and temporal data of the stable isotopes of water in groundwater, river water, and precipitation to determine seasonal bias in groundwater recharge and young streamflow in the Nelson River basin of west central Canada. Regionally in the northeastern US winter recharge ratios exceeds summer recharge ratios or the groundwater is biased by winter recharge (Jasechko et al., 2017). Cold-season recharge ratios were typically greater than warmer season recharge ratios and that precipitation which falls within the past 2 or 3 months make up about one-quarter of river discharge. Puntsag et al. (2016) examined a 43-year record of precipitation isotope values collected at the Hubbard Brook Experimental Forest in New Hampshire, US and noted a positive trend in the deuterium excess values and a negative trend in both the δ 18 O and δ 2 H. These trends were linked to increases in the interaction with air masses from the Atlantic Ocean and Arctic.
This study aims to examine the relationship between modern precipitation, surface water, and groundwater stable isotopes across Massachusetts (USA) via crowdsourcing to assess the isotopic variability of waters and relate variability to landscape and hydrological features. We also identify the spatial and temporal trends in seasonal surface and groundwater isotopes and 2-week weighted precipitation isotopes and distinguish hydrologic and hydrogeologic trends in surface and groundwater, respectively. Using these trends, we discuss and quantify the implications for either precipitation induced variability or variability due to open water systems, topographic differences, hydrogeologic and hydrologic processes in surface and groundwater stable isotopes. We envision that these datasets will allow us to address questions such as: (a) What is the nature of surface and groundwater interaction in the northeast US? (b) What are the potential impacts of climate change on stream flow generation, groundwater recharge, and groundwater storage? (c)How important are extreme precipitation events to groundwater storage? And (d) How does groundwater surface water interconnectivity change spatially across the landscape? Fundamentally, the use of stable water isotopes has become an inexpensive way to characterize the temporal and spatial variability of water bodies at a high resolution (Birkel et al., 2018) and can be used to inform management on both a local and regional scale.

BACKGROUND Stable Isotopes as Tracers
The interpretation and analysis of the composition of environmental water stable isotopes (precipitation, surface water, groundwater), δ 18 O and δ 2 H, are an important tool in examining the hydrologic processes on a global and regional scale (Dansgaard, 1964;Kendall and Coplen, 2001;Bowen, 2010;Puntsag et al., 2016). These analyses provide a better understanding for quantifying the spatially integrated effects of the water cycle and processes that occur in both the watershed and atmosphere (Bowen et al., 2011) as well as determining the relative amount of precipitation and groundwater in surface waters (Kendall and Coplen, 2001). Oxygen and hydrogen measurements of precipitation, surface water, and groundwater illustrate the effects of climate, topography, elevation, and various environmental parameters (Dansgaard, 1964;Lee et al., 2010;Welker, 2012;Evaristo et al., 2015;Akers et al., 2016).
Understanding the stable isotopes of water relies heavily on accurate and high precision measurements of δ 18 O and δ 2 H (Brand et al., 2009;Wassenaar et al., 2012). Isotope ratios are considered ideal tracers as they are part of the water molecule and can be easily sampled and preserved in groundwater, surface water, precipitation. Most importantly hydrogen and oxygen can preserve vital historical information (location, time, phase of precipitation) thus becoming a primary tool for hydrological, atmospheric, and meteorological studies (Reddy et al., 2006;Bowen et al., 2007;Timsic and Patterson, 2014). The stable isotopes of water are presented in the δ notation and represent the difference in heavy to light isotopes of water relative to the Vienna Standard Mean Ocean Water (VSMOW) (Craig, 1961;Sprenger et al., 2018). Although the delta notation is a dimensionless quantity the values are in per mil because of the low variation in the natural abundance of water stable isotopes (Kendall and Coplen, 2001;Sprenger et al., 2018 For the purposes of this paper, the term "enriched" will be used to describe water samples that have a high amount of heavy isotopes and "depleted" will be used to describe water samples that have a low amount of heavy isotopes. To determine the δ 18 O and δ 2 H of a water sample Equation (1) is used: where R is the abundance ratio of the heavy and light isotopes (e.g., 18 O/ 16 O and 2 H/ 1 H) and Rstandard is the VSMOW. In a dual isotope plot, δ 18 O -δ 2 H, the relationship between δ 18 O and δ 2 H is defined as the global meteoric water line (GMWL) (Craig, 1961) and is described by the following equation: isotopes. To determine the δ 18 O and δ 2 H of a water sample Equation (1) is used: This equation represents the relationship of δ 18 O and δ 2 H of surface waters globally and is an approximation of the mean world annual amount-weighted precipitation (Timsic and Patterson, 2014). This relationship is a result from Rayleigh processes, which is directly affected by temperature and pressure conditions during phase changes between liquid water and water vapor (Dansgaard, 1964). More recently, the stable isotopes of water are analyzed together with deuterium excess (d-excess), Equation (3), which was originally proposed by Dansgaard (1964).
d-excess is the y-intercept of the GMWL and is dependent on relative humidity, temperature, and kinetic isotope effects during evaporation (Kendall and Coplen, 2001). Because of this, d-excess values are sensitive to evaporative processes and can be used to measure the contribution of evaporated moisture and allow for additional assessments of environmental conditions during the time of vapor formation or rainout. High d-excess values indicate more evaporated moisture has been added and low values indicate samples fractionated by evaporation (Kendall and Coplen, 2001;Timsic and Patterson, 2014).
The line condition excess (lc-excess) is the deviation from the local meteoric line rather than the GMWL and is a good indicator of evaporative fractionation (Birkel et al., 2018).
It is defined by: Lc-excess = δ 2 H -a δ 18 O -b where a and b are the slope and y-intercept, respectively, of the GMWL.

Study Site and Region
Massachusetts is located in the northeastern portion of the United States and borders the Atlantic Ocean. It occupies 27,340 square kilometers where most of the state lies north of 42 • latitude In elevation, Massachusetts ranges from <152 to 1,063 m above sea level where the western portion is characterized as mountainous, the central portion as rolling hills and coast as flat land with marshes and small lakes and ponds (Boutt, 2017). Massachusetts lies in the prevailing westerlies, a region that is dominated by westerlies, generally eastward air movement and drier continental airflow (Weider and Boutt, 2010). Most of the precipitation events are sourced from colder regions: Arctic, Mid/North Atlantic and the Pacific. Massachusetts tends to see precipitation events that originate from the Arctic, Mid-Atlantic, North Atlantic, Pacific, Continental, and the Gulf of Mexico (Puntsag et al., 2016). Massachusetts receives ∼1,000 mm of rain annually and on average temperatures range from as high as 26 • C and low as −8 • C. This temperature variability is due to variance in topography. Because of this snowfall amount varies across the state making it difficult to determine a snowfall average. The National Climatic Division Center divided the state into three climatological divisions: Climate Zone 1, Climate Zone 2, and Climate Zone 3 (NCDC). Throughout this paper these zones will be shortened to CZ1, CZ2, and CZ3. Though, according to the Koppen-Geiger Classification, Massachusetts has four climatological divisions with the fourth climate zone encompassing the coast and Cape Cod. For this study, we have grouped together Climate Zone 4 with Climate Zone 3, thus dividing Massachusetts into only three Climate Zones.

METHODS
Most of the samples presented here were collected by volunteers representing different groups across the region. We provided and shipped volunteers with empty bottles, labels and standard sampling instructions. Samplers were instructed to fill bottles as full as possible and close the cap as tight as possible. Upon shipping back to us, we inspected all bottles and discard samples that were not more than 75% full (almost all bottles had no leaks). Repeat analyses monthly of stored samples show no evidence of evaporative enrichment. Inspection of isotopic results did not reveal any storage anomalies or sample irregularities. This was confirmed by replicate analysis and consistency of results from given locations or times of sampling. All isotope data presented here is reported in Supplementary Table 3.

Precipitation
Eleven sampling sites form the precipitation isoscape network presented here (Figure 1). Sampling for the precipitation network took place from March 2017 to March 2018 and is still ongoing. Each volunteer was supplied with thirty 30 mL high density polyethylene (HDPE) bottles, a small funnel, and a data collection sheet. Volunteers were asked to pour the contents of their rain gauge into a 1-L Nalgene bottle every day for a 2week time increment. At the end of the 2 weeks, the contents of the 1-L Nalgene bottle were poured into one 30 mL bottle which would result in a 2-week composite sample. During the winter months, if there was snow, volunteers would place the snow into a saucepan and dip the bottom of the pan into hot tap water. The melted snow can then be measured by pouring it from the saucepan into the inner cylinder of the rain gauge. Every 6 months, volunteers would send their precipitation samples to the University of Massachusetts-Amherst. Sample bottles were stored in plastic bags between the time of receipt and analysis.
Precipitation isotope analyses were screened and samples that had a deuterium excess <-10‰ were discarded as these may have been compromised by evaporation during storage/sampling. A total of 40 precipitation samples were removed from the analysis. After screening the precipitation isotope data we calculated a 2-week weighted averages following to remove any bias: Where P (2−weeks) is the 2-week precipitation amount as provided by the volunteers. We then determined the average annual with a 2-week weighted average following Where P i is the sum of the precipitation amount of the sampling site.

Surface Water Sampling Network and Data Sources
One thousand nine hundred and seventeen surface water samples from 556 surface water sites (Figure 1) across Massachusetts were analyzed for δ 18 O and δ 2 H. Surface water sites were selected based on their spatial location in order to accurately represent a surface water isoscape for Massachusetts. Surface water samples were collected from 2011 to 2018. Samples were collected from faculty and students from the University of Massachusetts-Amherst, watershed associations (Nashua Watershed, Charles River Water, Blackstone, Quabbin, and Chicopee River Watershed), and the acid rain monitoring (ARM) project. Volunteers were supplied with 1 or more 15 mL high density polyethylene (HDPE) bottles and a data collection sheet. On the day of sampling, volunteers were asked to thoroughly clean out the HDPE bottle with the to be collected water and then refill the bottle to the top to limit the amount of headspace and securely fasten the cap. Samples were returned to the University of Massachusetts where they were prepared for analysis within a few weeks upon arrival. Sample bottles were stored in plastic bags between the time of receipt and analysis.

Groundwater Sampling Network and Data Sources
One thousand four hundred and six groundwater samples from 409 groundwater sites across Massachusetts were analyzed for δ 18 O and δ 2 H (Figure 1 Each volunteer was supplied with 1 or more 15 mL high density polyethylene (HDPE) bottles and a data collection sheet. On the day of sampling, volunteers were asked to first run the tap water for a couple seconds, thoroughly clean out the bottle with the tap water and the fill the bottle to the top to limit the amount of headspace and securely fasten the cap. Samples were returned to the University of Massachusetts-Amherst where they were prepared for analysis within a few weeks upon arrival. Sample bottles were stored in plastic bags at room temperature between the time of receipt and analysis.

Analysis of Isotopic Composition of Water
Stable isotope compositions were measured by a wavelength scanned cavity ring-down spectrometry on un-acidified samples by a Picarro Cavity Ring Down Spectrometer (L2120-I) analyzer (Berden et al., 2000). Cavity ring-down spectroscopy is a direct absorbing technique which is conducted with eight pulse or continuous light sources which is significantly more sensitivity than the conventional absorption spectroscopy. The Picarro is equipped with a high precision vaporizer (A0211) and fitted with a CTC PAL auto-sampler. International reference standards (IAEA, 2000, Vienna, Austria) were used to calibrate the instrument to the VSMOW scale. To remove possible memory effect between samples, each sample was analyzed six times and the results of the first three injections were discarded. To further reduce memory effect, we adopted a modified version of a technique by Penna et al. (2012) where samples are grouped by water source and location. Three reference waters that isotopically bracket the sample values were run alternately with the water samples: Boulder, Colorado, Tallahassee, Florida and Amherst, Massachusetts, were used for a total of nine times each in every sample tray. The average values for these standards are −16.5, −2.6, and −7.5‰, respectively. These standards were calibrated with the Greenland Ice Sheet Precipitation (GISP), Standard Light Antarctic Precipitation (SLAP) and Vienna Standard Mean Ocean Water (VSMOW) from the IAEA (Kendall and Coplen, 2001). The results were calculated based on a rolling calibration so that each sample is determined by the three standards closest in time to that of the sample.

Precipitation
Three-hundred and ninety four (394) precipitation analyses are presented in Figure 2 on a dual isotope plot. δ 18 O values range from −18 to −1.3‰, δ 2 H values range from −132 to −8.0‰ and d-excess values range from −7.3 to 44‰. Average δ 18 O, δ 2 H, and d-excess values are −7.6 (±3), −47 (±24), and 12 (±6)‰, respectively. The precipitation data points form a flattened ellipse on the GMWL. Most of the data points fall on the GMWL with a few plotting off onto the evaporative enrichment line (Figure 2).
The State Meteoric Water Line (SMWL) generated from the unweighted stable isotope values of the 394 precipitation samples is δ 2 H = 7.7 * δ 18 O + 9.8. This regression line is calculated based on binning the 394 precipitation samples 0.5‰ in δ 18 O and then averaging the values. This is performed because the density of enriched samples are biased in the dataset, and the resulting line better represents the SMWL across the full range of isotopic compositions observed.
We compare our precipitation isotope analyses with two precipitation sites from Massachusetts collected by Vachon (2007). Weighted monthly averages from our network were calculated and compared to the weighted monthly averages of MA13-Boston, MA01-Cape Cod. Standard deviations between our weighted monthly averages and the weighted monthly averages of MA13 and MA01 ranged from 1.9 to 0.2‰.
Following the methods of von Freyberg et al. (2018), using the equation below, a sine-wave was fitted through the reported precipitation isotope values to better quantify the seasonal isotope cycles.
In Equations (6) and (7), A p is the amplitude for precipitation (‰), ϕ is the phase of the seasonal cycle, t is the time in years, is the frequency (yr −1 ), and k (‰) is a constant which describes the vertical offset of the isotope signal (von Freyberg et al., 2018). These equations allow us to quantify the amplitude of the seasonal isotope cycles and the coefficients, a, b, and c by using multiple linear regressions where a, b, and c are the amplitude, phase, and offset, respectively.
The amplitude A p is then determine by: Using these equations, we fit sine curves to the isotope data for all the precipitation isoscape samples, as well as each individual station (Figure 3). We found the equation of the sine fitted line to be: Cp(t) = −2.07 cos 2π × 1 365 x × −1.88 sin 2π × 1 365 x + −7.82. Precipitation varies approximately −5.4‰ seasonally. Amplitudes and offsets for each of the 11 precipitation sites were determined and compared with the calculated amplitudes and offsets of MA01 and MA13 (Vachon, 2007) to compare to our calculated values. The amplitude for each of our precipitation site ranges from 2.14 to 4.00 (±0.66) and the offset ranges from −6.12 to −9.39‰ (±0.99).
Additionally, seasonal local meteoric water lines were determined for winter and summer. Winter months includes October to March and summer months includes April to September. The LMWL for summer generated from the unweighted stable isotope values sampled during the summer months is δ 2 H = 6.6 * δ 18 O + 2.6. The LMWL for winter is δ 2 H = 7.4 * δ 18 O + 11. The slope, y-intercept, and d-excess of the winter LMWL is larger than the slope and y-intercept of the summer LMWL.
Median, minimum and maximum δ 18 O values were determined over a 12-month period and correlated to determine any statistical significance associated over this period, Figure 4. Over this 12-month interval, the δ 18 O values are variable where on average there is a −5.4‰ range. This same method was To assess the role of climatic variability on isotopic compositions, all precipitation samples were categorized into climate zones (Table 1). A trend toward more enriched precipitation values is observed moving from west to east across the study site (climate zones 1-3). Two-week weighted averages were determined for each of the 11 precipitation sampling sites, grouped into their respective climate zones and then related to their respective elevation (Figure 5). Elevations within the sampling sites range from 4 to 331 m. MA01 (Cape Cod) and MA13 (Boston) from Vachon (2007) were plotted alongside our sampling sites and categorized into CZ3. Sampling sites located at a higher elevations and in CZ1 and are more depleted than samples located at a lower elevation CZ3 which are more enriched.
Local meteoric water lines (LMWLs) were calculated for each of the climate zones. The climate zone 1 LMWL generated from the 31 samples located in that region is δ 2 H = 7.8 * δ 18 O + 13, climate zone 2 LMWL generated from the 186 samples located in that region is δ 2 H = 7.5 * δ 18 O + 9.3 and the climate zone 2 LMWL generated from the 178 samples located in that region is δ 2 H = 7.4 * δ 18 O + 8.1.
These maximum and minimum δ 18 O values are primarily driven by temperature and moisture source differences. National Oceanographic and Atmospheric Administration Air Resources Laboratory Hybrid Single-Particle Lagrangian Integrated Trajectory (HYSPLIT; Draxler and Rolph, 2003)    to 24.9‰, and lc-excess values range from −20.8 to 12.2‰. The dataset forms a flattened ellipse that lies slightly above and angled to the SMWL, Figure 6 illustrates the isotopic variability.  x + −8.48.
Using Equation (8), we found the amplitude of the sine curve to be 1.13‰, illustrated in Figure 8. Surface water experiences a seasonal range of 2.26‰ in δ 18 O.

Groundwater
The δ (Figure 9) with a higher percentage of points plotting off the SMWL and a few plotting off onto the evaporative enrichment line. Average δ 18 O, δ 2 H, and d-excess values for groundwater are 8.7 (±1), 53.3 (±6), and 16.3 (±4)‰, respectively. In order to better show the isotopic variability in groundwater, a histogram, which shows the frequency and distribution for both δ 18 O and δ 2 H were plotted and placed on appropriate axes with the dual isotope plot of groundwater analyses. At approximately δ 2 H −45‰ there is slight increase in the number of groundwater samples, this indicates where on the dual isotope plot groundwater samples start to fall off the SMWL and onto the evaporative enrichment line. The Groundwater Water Line (GWL), which is the unweighted linear regression generated from the 1,405 groundwater analyses is δ 2 H = 6.5 * δ 18 O + 2.9.
The slope of the GWL equation is 6.5 and the y-intercept of intercept is 2.9.
To visually represent seasonal variability of groundwater, the unweighted δ 18 O and δ 2 H values of the groundwater samples were compared with the time of sampling. Medians, minimum, maximum, 25th and 75th percentile, and 1% interquartile ranges were determined for each of the 12 months. Figure 10 illustrates little isotopic variability between the medians of each month and indicates relatively homogenous seasonal variations in the unweighted δ 18 O values. A sine-wave was fit through the reported groundwater isotope values to better describe the seasonal isotope cycles: Cp (t) = 0.05 cos 2π × 1 365 x × −0.04 sin 2π × 1 365 x + −8.79. Using Equation (8), we found the amplitude of the sine curve to be 0.07‰. Regionally groundwater experiences very little variability, with a seasonal range of 0.14‰.
Groundwater samples were grouped based on the aquifer they were sampled in: glacial till, outwash plain, bedrock, and glacial fluvial deposits. Out of the 1,405 groundwater samples the aquifer for 44 of the samples were unknown. Statistics of groundwater isotopic compositions are summarized in Table 1 and all data is plotted as a function of aquifer geologic type are shown in the dual isotope plot Figure 11. Samples located in outwash plains with a few samples from glacial fluvial aquifers plot off the SMWL. Outwash fluvial, and bedrock aquifers are isotopically very similar, while outwash plains are the most enriched.
Groundwater average δ 18 O values for each site were analyzed as a function of elevation to determine elevation effects. The elevation of groundwater samples ranges from 22 to 677 meters. Groundwater samples located at a lower elevation are more enriched than samples located at a higher elevation which are depleted. From 22 to 100 meters, δ 18 O values decrease from −3 to approximately −8‰. From 100 meters to 700 meters, δ 18 O values show some variability but are relatively similar at around −9‰.

Comparison of SMWL, SWL, and GWL and Isotopic Variability
Important differences in the composition of isotopes were observed between precipitation, surface water, and groundwater. The SMWL is δ 2 H = 7.7 * δ 18 O + 9.8, the SWL is δ 2 H = 5.7 * δ 18 O -4.2, and the GWL is δ 2 H = 6.5 * δ 18 O + 2.9, each line having its own unique slope and y-intercept (Table 1). Figure 12 illustrates the comparison between the SMWL, SWL and GWL and includes the 95th percent confidence interval for all three which was calculated using the curve fitting tool in MATLAB. The 95th percent confidence interval illustrates the minimum, maximum values and the range of variability for precipitation, surface water, and groundwater. These values and ranges allow us to compare and determine how variable the isotopic composition of δ 18 O and δ 2 H for precipitation, surface water and groundwater samples. The isotopic composition of precipitation across Massachusetts is more variable than groundwater and surface water. Both groundwater and surface water across Massachusetts experience some but smaller isotopic variability. Both the SWL and the GWL have shallower slopes than the SMWL and have an imbricate nature relative to the SMWL, Figure 12.
These differences reflect local precipitation processes but also reveal the correlation between precipitation, surface water and groundwater. Because the slope of the SWL and the GWL are both pulled toward the line of evaporative enrichment, it can be determined that the slope of both these lines are biased by evaporation in an open water system. The slope of the GWL is larger than the SWL but smaller than the SMWL. This is telling that groundwater is composed of some fraction of precipitation and surface. Although surface water is more variable than groundwater, the isotopic composition of groundwater have similar ranges, meteoric water line slope, and averages compared to surface water. Though in order to determine the fraction of surface water to precipitation in groundwater, a data comparison between modeled and observed analyses would have to be performed. As noted, the slope of the SWL is shallower than both the SMWL and the GWL. The isotopic composition of surface water is affected by both the isotopic composition of precipitation, the isotopic composition of groundwater, and open vs. closed water systems. As a result of this, the isotopic composition of surface water will be inherently weighted by precipitation and groundwater (Dutton et al., 2005) and it might be expected for the SWL to be more similar to the SMWL. From Kendall and Coplen (2001) low slopes of the LMWLs may be suggestive of post-rain evaporative enrichment which can be reflected in the surface water. The slope of the SWL is drawn downward and biased possibly by open water enrichment. In a study performed by Kendall and Coplen in 2001, LMWLs are likely to have slopes as low as 5 to 6 in arid regions where summer rains are primarily derived from the Gulf of Mexico and are the main source of recharge to surface water. Kendall and Coplen (2001) stated that many of the surface water samples, especially in arid zones, show LMWL slopes <6 are indicative of evaporation. Interestingly they found that most of the eastern sites in their study had slopes that range between 7 and 8 and have intercepts that ranged from 5 to 12‰. The slope of the SWL in Massachusetts is 5.7, which, as suggested by Kendall and Coplen (2001) indicates significant evaporation. As evidenced by the trend in isotope values in Figure 6, some surface water samples fall off and to the right of the meteoric water line (dexcess <10). Samples with these values are primarily from waters in the southeast part of the state with a preponderance of open surface water bodies.
Alternatively, one consideration that hasn't been discussed are anthropogenic affects such as the presence of dams or impoundments, which are present in some of the rivers and reservoirs included in this study. Dams naturally increase evaporation, thus rivers or reservoirs associated with these dams will produce an evaporative signal and cause the δ 18 O values of surface waters associated with these dams to be high. This in effect will cause the recycling of surface water that already have high δ 18 O values to move through the hydrologic cycle where further enrichment of δ 18 O values may take place. Another consideration is the order of the rivers and/or streams and the size of the catchment of the rivers. When we compare the slope of the GWL to the SMWL and the SWL, we find that it lies between both the SMWL and the SWL. With the exception of <5% of the groundwater samples, most of the samples lie above the SMWL indicating high d-excess values. Those samples that have δ 2 H and δ 18 O values that lie below the SMWL and onto the evaporative enrichment line are indicative of samples that have gone through isotopic fractionation. Most of the groundwater samples that have an evaporative signal are located in CZ 3 where outwash plain aquifers are located and are only found in the southeast portion of Massachusetts. Outwash plains in the southeast portion of Massachusetts have a high amount of surface water features such as kettles lakes and ponds where ponds are essentially the water table. This could cause a high amount of groundwater surface water variability. As demonstrated in Figure 11, outwash plains are the most enriched and also the most variable. This portion of Massachusetts was located near the terminus of an ice sheet and consists of sandy terminal moraines. Because of the permeable nature of the aquifer, an estimated 45% of the average annual precipitation percolates into the soil and becomes groundwater recharge and an estimated 55% of the precipitation is evaporated (Olcott, 1995). Figure 13A illustrates the distribution and variability of δ 18 O values for weighted precipitation, surface water, and groundwater. The median for groundwater, surface water, and precipitation are −8.78, −8.16, and −6.88‰, respectively. The positive ends of the surface water and groundwater plot are influenced by more enriched precipitation. This offset reveals that the isotopic signal of precipitation becomes filtered due to recharge infiltration, evaporation, and evapotranspiration and is reflected in the surface and groundwater. The depleted values in both the surface water and the groundwater are suggested to be primarily due to snowmelt and winter precipitation events.
Lc-excess was calculated for both groundwater and surface water to define the number of groundwater and surface water samples that have gone through evaporative enrichment. Negative lc-excess values indicate evaporative enrichment and positive lc-excess values indicate differences in moisture source (Birkel et al., 2018). From Figure 13B it can be shown that approximately one quarter and one eighth of surface water and groundwater samples, respectively, have lc-excess values less than zero. Though it is no surprise that there are more surface water samples with higher lc-excess values than groundwater samples because surface water samples experience more isotopic fractionation than groundwater. We found that the SMWL for precipitation to be almost identical to the global meteoric water line. The global meteoric water line (GMWL) is defined as δ 2 H = 8 * δ 18 O + 10 and the SMWL is δ 2 H = 7.7 * δ 18 O + 9.8. When compared to the GMWL, the slope and y-intercept is less than the GMWL. Differences in slopes could be to in-cloud processes such as sub-cloud evaporation, evaporation during rainfall, or intensity of rainfall (Ren et al., 2016).

Seasonal and Hydrologic and Hydrogeologic Effects
Precipitation Precipitation experiences a seasonal isotopic enrichment and depletion in the heavy isotopes. During the summer and spring months storms are primarily sourced from areas where temperatures are high and the isotopic composition of water vapor are more enriched. The moisture that becomes precipitation is usually sourced from surface water at mid to low latitude oceans where elevation is low and temperatures are high and evaporation occurs from the continents (Welker, 2012). During the winter months snow is the primary mode of precipitation. Massachusetts sees Nor'easters, which are storms that are sourced from the Mid-Atlantic and are formed due to the sharp contrast in temperature between the Gulf stream current and the cold air masses from Canada. Due to the formation of these storms, the storms will produce isotopically depleted precipitation. Typically, precipitation events during winter months are sourced from surface water at high latitudes and high elevations where temperatures are cold. For the most part precipitation in the Northeast is advected from the Gulf of Mexico, Mid-Atlantic, Pacific, Arctic, North Atlantic, and Continental (Welker, 2000).
Many researchers have determined and documented significant differences between the isotopic compositions of winter vs. summer precipitation (Birkel et al., 2018, von Freyberg et al., 2018. Birkel et al. (2018) determined the LMWLs for fall, winter, spring, and summer for 22 catchments in Scotland and found that between all the seasons winter had a larger slope than summer. Winter and summer meteoric water line were determined for the 394 precipitation samples in our precipitation isoscape network. Similar to the findings of Birkel et al. (2018), the slope of the winter meteoric water line is larger than the slope of the summer meteoric water line.
These differences are primarily due to seasonal differences in stormtracks, but can also be due to seasonal changes at the precipitation site (Kendall and Coplen, 2001) and reflect changes in source and temperature during a hydrologic year. In Massachusetts, there is a seasonal temperature discrepancy between summer and winter months. Average monthly temperatures were collected from the NOAA National Centers for Environmental information, Climate at a Glance: Global Time Series. The average temperature for summer, which is defined as April to September (Jasechko et al., 2014) is 16 • C and average winter, which is defined as October to March (Jasechko et al., 2014) temperature is 3 • C (NCDC). It is documented that air temperature and average 2-week weighted monthly δ 18 O values have a positive correlation. As temperature increases precipitation becomes isotopically enriched and vice versa. At higher temperatures, the variability is less compared to lower FIGURE 12 | 95% confidence interval comparison between the SMWL, SWL and GWL. The arrows indicate low to highest values within the 95% confidence interval polygon.
temperatures where there is a larger amount of variability. This is primarily due to seasonal moisture source differences. HYSPLIT models shows that in the first week of March 2018, all sampling sites experienced precipitation events that originated from a continental source, had a more northerly route and experienced cooler temperatures and more moisture recycling due to the fact that the air mass had a continental pathway.

Surface Water
A seasonal trend was documented for our entire surface water dataset. During the summer months (April to September) surface waters experienced isotopic enrichment and during the winter months (October to March) surface waters experienced an isotopic depletion, Figure 7. April experienced the lowest δ 18 O values while September experienced the highest δ 18 O values. Interestingly, as mentioned above summer is defined by the months of April to September and winter is defined by the months of October to March; though it is the month of April that experiences the most isotopic depletion. One possible explanation for this difference is differential recharge or delayed snowmelt (von Freyberg et al., 2018). The average temperature of March in Massachusetts is ∼0 • C while the average temperature of April is ∼10 • C. This temperature difference could cause snow to melt much later in the winter season and into early spring. Because of the late snowmelt, surface water would experience recharge that has lower δ 18 O values. Dettinger and Diaz (2000) found that in the eastern United States the amount of precipitation entering streams is higher during the winter months than the summer months.
Modifications in isotope composition is impacted by evaporative enrichment during summer months, the amount of enrichment due to surface water type (stream vs. pond/lake), and the source of precipitation all affect the slope the SWL. δ 18 O comparisons for each sampling month show a large number of surface water samples experience enrichment in the heavy isotopes and plot above the 1% interquartile range. One possible explanation for these outliers are that these samples experienced open water enrichment isotopic fractionation causing sample to become isotopically enriched. Alternatively, these samples could be the result of additions of isotopically enriched precipitation. Samples categorized based on the type of surface water show that the δ 18 O values of each type lakes, reservoirs, and ponds were more variable (and show more d-excess and enrichment off of the SMWL) than rivers, streams, and springs.
In order to investigate whether precipitation or open water enrichment are the primary effects that cause surface water isotopic variability we used the amplitudes calculated from the sine-wave fit lines in both precipitation and surface water to determine the fraction of young water to precipitation. As defined by von Freyberg et al. (2018) the young water fraction is the proportion of catchment discharge that is younger than 2-3 months and can be approximated from the amplitudes of seasonal cycles of stable water isotopes in precipitation and surface water. Young water fractions have been inversely correlated with water table depth and topographic gradients (Jasechko et al., 2014) and have become a useful value for catchment inter-comparison as well as catchment characteristics as it can be calculated from limited data (von Freyberg et al., 2018). By looking at the seasonal isotopic composition of precipitation one can determine the damping and phase shift of the seasonal cycle as it gets propagated through the catchment which can be used to determine the timing of catchments (von Freyberg et al., 2018). Higher stream flows will have larger young water fractions as an increase in stream discharge due to rain events will contain more recent precipitation than base flow. Lower stream flows will have smaller young water fractions and will contain more base flow.
From our data we find the amplitude of the mean δ 18 O of precipitation is 2.7‰ and the mean amplitude of δ 18 O of surface water to be 1.13‰. To calculate the fraction of young water to precipitation we used the following equation (von Freyberg et al., 2018) below where Ap and As are the amplitude of precipitation and surface water respectively: Fraction of young water = As/Ap. From our date we calculate the average fraction of young stream water to be 0.4. This value indicates that ∼40% of total discharge is composed of "young waters" or water that is <2-3 months in age.
Surface water follows the same seasonal enrichment and depletion as precipitation but dampened, with a smaller amplitude (Figure 8). The phase shift between precipitation and surface water is very small, almost non-existent. This is likely due to a delayed release of depleted winter precipitation from snowpack (von Freyberg et al., 2018). A smaller damping also implies a smaller fraction of young water in the stream flow, a smaller proportion of modern water to base flow (von Freyberg et al., 2018). The lack of phase shift between the seasonal cycles of surface water and precipitation along with the calculations to determine the young water fraction, implies that surface water has a higher proportion of modern water to base flow.
Residual corrected model values compiled from Bowen et al. (2011)  Many of the data points lie within the seasonal variability range, suggesting a strong role of evaporation enrichment. P:PET induced seasonal variability is a primary control on surface water variability. Data points that plot outside of the seasonal range is indicative of open system enrichment (i.e., evaporative enrichment). These results suggests that 90% of surface water variability is due to P:PET induced seasonality and 10% is due to open system induced variability.

Groundwater
Global seasonal groundwater recharge (Jasechko et al., 2014) calculated by using the seasonal differences in groundwater recharge ratios between the summer and winter seasons are helpful for assessing water resource dynamics. Based on the calculated groundwater recharge ratio (Jasechko et al., 2014), recharge in eastern Massachusetts have winter time groundwater recharge ratios that are consistently higher than summertime recharge ratios. Temperate regions, which includes eastern Massachusetts, that a large percent of winter precipitation infiltrates and recharges aquifers relative to summer precipitation. Our data presented here (based on many samples) suggests larger magnitude of bias toward winter time recharge. For example the mean δ 18 O of precipitation we report here is −7.6‰ compared to groundwater mean δ 18 O of −8.7‰ (Table 1). Jasechko et al. (2014) reports similar mean isotopic composition of precipitation of −7.5‰ with a groundwater composition of −7.9‰. This bias is not necessarily even across the 3 climate zones. Strong variations in the isotopic composition of the more topographically varied climate zone 1 may explain the fact that the mean δ 18 O of precipitation in climate zone 1 being −9.3‰ compared to −9.1‰ for groundwater may suggest some landscape control on compositions and inferred differences. Additionally, the dataset does suggest that latesummer growing season high intensity precipitation events may be an important contributor to groundwater recharge (e.g., Boutt et al., 2019).
The four types of geologic deposits in the region (glacial till, glacial fluvial, bedrock, and outwash plains) have unique characteristics, properties and isotopic compositions. Outwash plains have the highest average δ 18 O value, glacial till, glacial fluvial, and bedrock have similar average δ 18 O values and are the most depleted. Most of the enriched samples plot higher than the 1% interquartile range, suggesting these samples are most likely groundwater that have experienced some sort of fractionation. Outliers that plot below the 1% interquartile range δ 18 O value could be representative of the isotopic composition of direct recharge via snow melt. Groundwater located in outwash plains are more variable than bedrock, glacial fluvial, and glacial tills. In the regions bedrock aquifers, which are dominated by fractured crystalline rock, they do not recharge quickly as the movement of water is controlled by the existence of secondary openings in the aquifer, which can create long residences time flow paths (Weider and Boutt, 2010). Groundwater in bedrock aquifers have FIGURE 14 | Average surface water, groundwater δ 18 O cross plot of samples binned by HUC10 watershed location. Surface waters within the HUC10s tend to be more enriched compared to groundwaters. a long residence time because of its low permeability. Till aquifers are located at higher elevations and primarily in recharge areas (Weider and Boutt, 2010). They tend to be thin and consist of unconsolidated material and drain and fill seasonally due to the permeability and porosity of the aquifer material. Stratified glacial fluvial aquifers are found in valleys and areas near streams. They tend to have a larger storage capacity than till aquifers (Weider and Boutt, 2010). These hydrogeologic properties as well as their water storage capacity are the determinant factors of groundwater and surface water interaction and may be an implication as to the lack of seasonal variability seen in the groundwater δ 18 O values.
Unlike precipitation and surface water, groundwater does not experience a statistically-observable seasonal isotopic enrichment and depletion due to what is known as a reservoir effect which indicates a well-mixed reservoir of groundwater. The young water fraction for groundwater was also determined and calculated to be ∼0.02‰ which indicates that 2.60% of groundwater is <2-3 months in age. This lack of seasonal variability can also be due to a small seasonal input which can be further dampened by groundwater residence time being higher than the volume of the groundwater storage. δ 18 O values for each sampling month were compared to one another. In general, the groundwater δ 18 O values for each sampling month does not show any significant enrichment or depletion during the summer and winter months, respectively, but rather the δ 18 O values are fairly constant. This lack of seasonal variability may suggest that the groundwater in Massachusetts is being dominated by processes occurring in the vadose zone, or local variability due to hydrogeologic properties.
Spatial variability in isotopic compositions may be due to a high surface water to groundwater interaction which may be attributed to aquifer characteristics. To investigate this, we determined mean groundwater and surface water sample sites were grouped by HUC10. Surface water average δ 18 O values are compared to groundwater δ 18 O values in Figure 14. A one to one line is plotted to determine the percent of points that either plot above or below the line. If the isotopic composition of groundwater and surface water were the same, all points would plot on the one to one. Most (82.5%) HUC10 basins fall to the left of the one to one line, which indicates that a majority of the surface water samples are more enriched than the groundwater samples. This further suggest an important role of evaporative  (1)

(3)
All values reported are means. enrichment or seasonal bias in surfacewater compared to shallow groundwater.

CONCLUSIONS
Isoscapes and isotopic studies are an important tool to better understand hydroclimatic studies and their effects on water resource and water management. The use of stable isotopes has become a low cost, effective tracer mechanism to gain more knowledge on precipitation dynamics, groundwater recharge mechanisms, paleoclimate, and water resource management (Birkel et al., 2018). A dataset of 1,707 surface water samples, 1,405 groundwater samples, and 558 precipitation samples across Massachusetts has been analyzed in terms of seasonal, temporal, spatial, and environmental variability with the aim of determining and explaining the isotopic signature and variability of precipitation, surface water and groundwater for the state of Massachusetts. A state meteoric water line, surface water line, and groundwater line have been calculated from unweighted precipitation, surface, and groundwater analyses. The SMWL is δ 2 H = 7.7 δ 18 O + 9.8. The SWL is δ 2 H = 5.7 δ 18 O + 4.2, and the GWL is δ 2 H = 6.6 δ 18 O + 4.0, Table 2. SMWL, SWL, and GWL reveal significant differences in slopes and yintercept, slopes range from 0.9 to 2.0‰. Mean δ 18 O and δ 2 H values of precipitation vary as a function of moisture source, seasonality, elevation and topographic location. Mean δ 18 O and δ 2 H values of surface water vary as a function of open vs. closed water systems, topographic location, seasonality, and precipitation induced seasonality. Mean δ 18 O and δ 2 H values of groundwater vary as a function of hydrogeologic characteristics and topographic location.
Spatially, δ 18 O and δ 2 H values for precipitation, surface water, and groundwater illustrate an isotopic separation along an east-west topographic gradient where isotopes are enriched in CZ3, where there is little change in elevation, and depleted in CZ1, where topography varies more. This east to west depletion is primarily due to differences in elevation and relative distance from a moisture source. Median δ 18 O values from 2017 and surface water samples 2018 for CZ1, CZ2, and CZ3 were correlated with the median δ 18 O values of precipitation. Precipitation median δ 18 O values plot consistently higher than both groundwater and surface water datasets. These results imply that the isotopic signature of precipitation gets filtered due to recharge infiltration, evaporation, and evapotranspiration and is reflected in both the surface water and groundwater.
Estimates of young water fractions for surface waters reveal that roughly 40% of total discharge is composed of water that is <2-3 months and that 2.5% of groundwater consists of water that is <2-3 months in age. Correlation between river water data that incorporates monthly precipitation isotopic variation illustrate that roughly 90% of our surface water variability is due to precipitation-evapotranspiration induced seasonality and 9% is due to open water losses. In this paper, we have developed a basic characterization and understanding of the isotopic signatures and variability of precipitation, surface water, and groundwater. We have related the interaction between surface water and groundwater and surface water to precipitation and have determined that a large portion of surface water variability is due to precipitation induced seasonality and groundwater variability reflects the dampening of this seasonality. The good agreement of the spatial patterns in precipitation, surface water, and groundwater will be useful for analyzing the effects of changes in moisture sources and timing and how these changes are propagated through the hydrological cycle in Massachusetts. Furthermore, by relating the surface water to precipitation and surface waters sensitivity to precipitation isotopic variability it will be useful in analyzing the changes in moisture source and how it is reflected in the stable isotope of surface water. This work has the potential to answer many more questions about the hydrologic cycle of Massachusetts and how factors such as catchment size, recharge processes in the vadose zone, and climate change can affect the isotopic topic signature of waters throughout the state of Massachusetts. Further work into studying these factors will improve our basic understanding of the hydrological behavior in a humid post-glacial environment.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Materials, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
DB initiated the study, designed and initiated the sample collection, proposed the research questions, secured funding, and led final compilation of manuscript. AC was primarily responsible for data analysis, figure creation, and writing of manuscript. All authors contributed to the article and approved the submitted version.