Footprint of recycled water subsidies downwind of Lake Michigan

. Continental evaporation is a significant and dynamic flux within the atmospheric water budget, but few methods provide robust observational constraints on the large-scale hydroclimatological and hydroecological impacts of this ‘ recycled-water ’ flux. We demonstrate a geospatial analysis that provides such information, using stable isotope data to map the distribution of recycled water in shallow aquifers downwind from Lake Michigan. The d 2 H and d 18 O values of groundwater in the study region decrease from south to north, as expected based on meridional gradients in climate and precipitation isotope ratios. In contrast, deuterium excess ( d ¼ d 2 H (cid:2) 8 3 d 18 O) values exhibit a significant zonal gradient and finer-scale spatially patterned variation. Local d maxima occur in the northwest and southwest corners of the Lower Peninsula of Michigan, where ‘ lake-effect ’ precipitation events are abundant. We apply a published model that describes the effect of recycling from lakes on atmospheric vapor d values to estimate that up to 32 % of recharge into individual aquifers may be derived from recycled Lake Michigan water. Applying the model to geostatistical surfaces representing mean d values, we estimate that between 10 % and 18 % of the vapor evaporated from Lake Michigan is re-precipitated within downwind areas of the Lake Michigan drainage basin. Our approach provides previously unavailable observational constraints on regional land-atmosphere water fluxes in the Great Lakes Basin and elucidates patterns in recycled-water fluxes that may influence the biogeography of the region. As new instruments and networks facilitate enhanced spatial monitoring of environmental water isotopes, similar analyses can be widely applied to calibrate and validate water cycle models and improve projections of regional hydroecological change involving the coupled lake-atmosphere-land system.


INTRODUCTION
Land-atmosphere fluxes of water and latent heat are significant and spatially heterogeneous components of the continental climate system that are difficult to measure over large spatial scales (Huntington 2006). These processes can be particularly profound in the vicinity of large lake systems, where they exert direct and indirect influence on the climate and ecology of surrounding land areas. Many such effects have been recognized and studied within the Great Lakes region of North America, including but not limited to lake effects on surface temperature, wind speed, precipitation amount and seasonal distribution, soil development, and distribution of forest types (Scott and Huff 1996, Schaetzl 2002, Henne et al. 2007. Robust characterization of land-atmosphere water fluxes in general, and lake effects in particular, is critical given the potential for significant changes, including amplification and non-linear feedbacks, in response to anthropogenic climatic and environmental alteration. Among the changes recognized within the Great Lakes region in recent decades are shifts in seasonal patterns of lake level change (Lenters 2001), timing of ice melt and surface water temperatures (Austin and Colman 2007), lake effect snowfall amounts (Braham andDungey 1984, Burnett et al. 2003), and lake water evaporation rates (Assel et al. 2004, Hanrahan et al. 2010. Although historical meteorological records provide a basis for the study of many lake effects on regional climate, one aspect of lake/land connectivity that is not well understood is the direct effect of large lakes on regional precipitation through recycling of water (i.e., lake-surface evaporation and subsequent re-precipitation). Classic studies of 'lake effect' precipitation have analyzed spatial patterns of precipitation rates downwind of large lakes, but these patterns integrate effects associated with both moisture and energy transfers between the lake surface and atmosphere (Changnon and Jones 1972, Braham and Dungey 1984, Scott and Huff 1996. In contrast, studies using the stable isotopes of H and O as tracers can deconvolve these effects because addition of evaporated lake water produces a predictable increase in the deuterium excess values (d ¼ d 2 H À 8 3 d 18 O) of atmospheric water vapor and subsequent precipitation (Gat 1996).
Previous applications of stable isotope tracers to study lake water recycling in the Great Lakes region have suggested that up to 15% of summertime atmospheric vapor and precipitation sampled downwind of the lakes is derived from lake water recycling (Gat et al. 1994, Machavaram andKrishnamurthy 1995). Traditionally, monitoring of meteoric water at high spatial resolution has been logistically challenging, and it has not been possible to evaluate spatial variation in recycled water contributions to terrestrial environments, estimate the magnitude of this flux over hydrologically significant areas, or assess the relevance of recycling to the water budget of the Great Lakes Basin. In this study we use a spatially distributed network of groundwater sampling sites in order to quantify the spatial footprint and regional magnitude of recycled water fluxes. Our work demonstrates a novel GIS-based approach to the quantification of recycled water fluxes, which can be extended using other water isotope monitoring data to provide previously unavailable observational constraints on regional recycling dynamics.

METHODS
We collected water samples from 108 private wells during the period February 2008 to August 2009 by running cold tap water for 30 seconds prior to filling, capping and sealing a 1 dram vial with poly-lined seal. Target wells were selected from a database maintained by the Michigan Department of Environmental Quality groundwater mapping project. Selection criteria targeted recently recharged water, including wells completed in unconfined aquifers with a screen depth of less than 100 feet and situated on coarsetextured glacial outwash deposits. Samples were refrigerated at 48C prior to isotopic analysis, and were analyzed at the Purdue Stable Isotope Laboratory (http://www.purdue.edu/eas/psi/) within approximately 3 months of collection. Isotopic analysis was conducted in triplicate using a Thermochemical Elemental Analyzer coupled to a Delta V mass spectrometer (Ther-moFinnigan) and previously published methods (Gehre et al. 2004). All data were normalized to the Vienna standard mean ocean water-standard light Antarctic precipitation (VSMOW-VSLAP) scale through analysis of 2 laboratory reference waters and are reported as d values relative to the VSMOW standard (Coplen 1996). The average analytical precision was 0.7% for d 2 H and 0.2% for d 18 O (1 r) based on the repeated analysis of a third reference water, giving an average precision for d values of 1.1% (assuming that errors in the H and O isotope ratio measurements of a given sample were independently distributed).
Data were analyzed relative to gridded precipitation data  obtained from the Parameter-elevation Regressions on Independent Slopes Model (PRISM) project (Daly et al. 1994) and interpolated depth of freshly-fallen snow (Daly and Taylor 2000). Climate variable values were extracted at groundwater cell sites using Hawth's Tools in ArcGIS 9.3 and analyzed relative to the isotopic data in Microsoft Excel 2007 using the Analysis ToolPak.
We estimated the fractional contribution recycled Great Lakes water to the groundwater samples using the model of Gat et al. (1994), which describes the d-separation between atmospheric moisture traversing a lake (d a ) and recycled moisture added from lake-surface evaporation (d E ) as: where d w is the deuterium excess value of lake surface water, h is the relative humidity of the boundary layer over the lake, the constant 107 reflects the effect of kinetic fractionation on d values, and h describes the perturbation of atmospheric humidity by lake evaporation. Using this expression, the fractional contribution of recycled water ( f r ) to a water sample taken down-wind of the lake system and having a deuterium excess value d s is (assuming no subsequent non-equilibrium fractionation): Evaluation of the spatial structure of the groundwater data was conducted using the Moran's I test (Moran 1950). Computation of this test statistic and subsequent interpolation of isotopic values throughout the study region was conducted using a new automated geostatistical code written by the authors for use in the IsoMAP web-GIS modeling portal (http:// isomap.org). This code was used because it supports geostatistical interpolation with ancillary independent predictor variables and because it allowed efficient iterative model-fitting for model development and uncertainty analysis. Source code and comprehensive documentation are available in Ecological Archives, but briefly we model the isotopic values using the geostatistical model: where Y is the vector of isotopic values at observed locations, X is the matrix of independent variable values, b is the vector of parameter coefficients, and e is the vector of error terms. In order to account for spatial dependence in e (as inferred from the Moran's I test) we describe spatial correlation in e with the Matérn correlation function (Matérn 1986). The Matérn function has three parameters k ¼ (k 1 , k 2 , k 3 ): a strength parameter k 1 that controls the magnitude of spatial dependence, so that Eq. 3 reduces to a linear regression model if k 1 ¼ 0; a range parameter k 2 that controls the rate of decay of spatial dependence; and a smoothness parameter k 3 that controls the behavior of the smoothness of the correlation function. The code numerically optimizes the likelihood function of Eq. 3 and derives the estimates of the model parameter coefficients (b) and Matérn correlation function parameters (k). These are used to calculate the predicted value at location s 0 using the universal kriging predictor (Stein and Corsten 1991): where x 0 0 is a vector of independent variable values at s 0 and c 0 0 and Rk À1 are the vector and matrix of modeled correlation values for Y * (s 0 ) relative to all measurement sites and for e, respectively. Input includes a table of coordinates (latitude and longitude), observed isotopic values, and independent variable values at measurement sites, and an equivalent table (excluding the isotopic values) for all sites at which Y * is desired. Output includes summaries of all model parameters, including standard deviations and significance (t test), as well as tables of predicted isotopic values and the standard deviation of these predictions.
For this work we used the statistical code both to generate prediction surfaces using our entire dataset and to conduct leave-one-out crossvalidation. For cross-validation, the dataset was iteratively re-sampled and the code used to optimize the model on the partial dataset and predict values both at the excluded observation site and at grid cells across our spatial study domain. Predictions were made on a 0.05 3 0.058 grid using a model parameterized in terms of site latitude or longitude, and were clipped to land areas within 50 km of a measurement site in ArcGIS 9.3.
Because we are interested in linking the groundwater samples to atmospheric processes, we evaluated the observed geographic distribution of values in the context of documented regional precipitation isotope ratio patterns. In most non-arid climates, the isotopic composition of shallow unconfined aquifer water is similar to that of the amount-weighted annual average precipitation in the recharge area (Gat 1996, Smith et al. 2002, but potential exists for postprecipitation modification of these values, particularly through evaporation. Evaporation of meteoric water produces an increase in the isotope ratios and a decrease in the deuterium excess value of the residual water (Gat 1996, Clark andFritz 1997). Most samples (91 of 108) have d values higher than þ10% and lie above the global meteoric water line (Fig. 1A). Ten samples have d values that are much lower (less than þ8.0%), plotting well below the GMWL and below values characteristic for regional precipi-tation (Gat et al. 1994, Bowen andRevenaugh 2003), and likely reflect evaporative effects. We also observe strong correlation between the hydrogen and oxygen isotopic compositions of the water samples and site latitude (Fig. 1B), as is typical for mid-latitude continental precipitation (Rozanski et al. 1993). We identify nine 'outlier' samples, however, which have d 2 H values that are higher than the trend defined by the other data by 10% or more (Fig. 2). Of these nine samples, 8 also have d values lower than þ8.0%. Based on these two criteria, we exclude data from a total of 11 samples, likely affected by appreciable evaporative fractionation, from further analysis.
Isotopic data within the screened dataset are significantly correlated with a number of geographic and climatic parameters ( Fig. 1, Table 1). Hydrogen and oxygen isotope ratios are strongly and significantly correlated with latitude and longitude, as well as with distance from the Lake Michigan coast (shortest straight-line distance), mean annual temperature (MAT), precipitation (MAP), and snowfall (MAPS; each of which is autocorrelated with the geographic variables). Relationships between sample isotopic composition and latitude and temperature are particularly strong, and the slope of the relationship between groundwater d 2 H values and MAT (À6.0%/8C) is similar to the 'global' mean annual precipitation d 2 H/MAT relationship for midlatitude sites (À5.6%/8C, (Dansgaard 1964)). Together, these observations suggest that the screened groundwater dataset likely represents a robust proxy for variation in mean annual precipitation isotope ratios across the study region.
Deuterium excess values, in contrast, are correlated with longitude, distance from the Lake Michigan coast, and the two precipitation variables but not latitude or temperature (Table 1, Fig. 1C), suggesting that climatological controls on d and d 2 H are decoupled. Values decrease from an average of þ15.8% for sites within 10 km of the Lake Michigan coast to þ13.1% for sites more than 100 km from the coast. Groundwater d values were weakly correlated with snow amounts. High deuterium excess values are characteristic of some solid-phase condensate due to kinetic effects accompanying the rapid growth of ice crystals under supersaturated v www.esajournals.org v www.esajournals.org conditions (Jouzel and Merlivat 1984). The existence of correlation between d and snowfall amount across the study region strongly suggests that spatial variation in kinetic effects associated with solid-phase precipitation could contribute to the observed groundwater d patterns. However, the coefficient of determination for this relationship suggests that variation in snowfall amount explains only 7% of the observed variation in d values. The lack of correlation between the observed d values and temperature, and the relatively weak correlation between d and MAP, indicate that patterns of d values in our data set are not strongly linked to variation in variables related to potential evaporation rates.
Although we cannot completely rule out the influence of post-precipitation evaporation as a control on the groundwater d values in our screened dataset, this observation suggests that it is unlikely that evaporation is the dominant control on the large-scale deuterium excess patterns observed here.

Recycling model
The high d values of Michigan groundwater samples are consistent with the interpretation that these samples include a significant amount of water recycled from the Great Lakes, as previously documented for precipitation, stream water, and vapor samples in the region (Gat et al.   (Table 2), and þ1.8 for Lake Michigan surface water based on collections of Lake Michigansupplied tap water (Bowen et al. 2007), supported by direct measurements of lake-water d 18 O (Dufour et al. 2005). Gat et al. (1994) estimated the values of h to be 0.88 for summer conditions over Lake Superior, but emphasized that this was the most poorly constrained parameter in their model. We used the formulation of Gat et al. (1994, Gat 1995 to re-calculate monthly values of h for Lake Michigan based on climatological normal surface-air temperature and dew point data from Sheboygan, WI and U.S. National Data Buoy Center Station 45007 (http://www.ndbc.noaa. gov/). Data for individual months (value were only available for March-December) give a wide range of values from 0.16 to 1.14 (effectively ¼ 1.00), but center on an average value of 0.61 with no strong seasonal trend (Fig. 3). This relatively low value suggests a somewhat larger perturbation of atmospheric humidity over Lake Michigan than estimated by Gat et al. (1994) for Lake Superior. This may in part reflect higher evaporation rates over Lake Michigan, but may also be influenced by the data used: for example, the stations compared here for Lake Michigan are separated by a degree of latitude and may at times be affected by different airmasses. To account for this uncertainty, we calculate f r for values of h ranging from 0.61 to 0.88. Although a comprehensive analysis of uncertainty in other model parameters is beyond the scope of this manuscript, we note that in general these are better constrained than h, and available data for other potential water sources to the regional atmosphere (e.g., Lake Superior) indicate similar values for most parameters (e.g., surface water d values for Lake Superior ' þ3.0% (Mae and  Longstaffe 2006)). Overall, this suggests that uncertainty associated with our estimates for these parameters is minor relative to that associated with h.
Our calculated values for f r at individual groundwater sites range from 0 to 32% (average ¼ 14%) and from 0 to 18% (average ¼ 8%) using the lower and higher values of h, respectively. These estimates encompass previously published estimates of f r for atmospheric vapor and precipitation in the region (Gat et al. 1994, Machavaram andKrishnamurthy 1995), and suggest that recycled lake water comprises a significant component of subsurface freshwater resources in Michigan's Lower Peninsula. Recycled water fractions for the measurement sites vary spatially in parallel with the measured values of d. The average f r estimated for sites west of À86.08 longitude is approximately 60% higher than that of sites in the eastern part of the study area (longitude . À84.58; t test, P , 0.01), implying that the recycled water flux is of greater relative significance to hydro-and ecosystems in areas closer to the Lake.

Spatial analysis
Groundwater H and O isotope ratios measured in this study show significant spatial autocorrelation after detrending for the observed linear relationship with site latitude (Moran's I ¼ 0.14, P , 0.001). We used the IsoMAP geostatistical code to produce an interpolated map of d 2 H values throughout the study domain using site latitude as a predictor variable (Fig. 4A). The mean absolute error (MAE) and root mean squared error (RMSE) values for these predictions, calculated from the cross-validation results, are 2.6% and 3.5%, respectively. These values represent approximately 6% and 9% of the total range of observed values, and are much smaller than the uncertainties typically associated with global water d 2 H models (Bowen and Revenaugh 2003). This result is consistent with other studies that have shown that region-specific models can provide enhanced water isotope prediction accuracy compared to global analyses (Hobson et al. 2004, Dutton et al. 2005, Lykoudis and Argiriou 2007, Lechler and Niemi 2011. Our prediction map depicts strong, dominantly meridional,  Dutton et al. 2005). Additional structure is apparent in some parts of the map. In particular, the d 2 H values at most sites in close proximity to the Lake Michigan coastline tend to be several % lower than those at sites further to the east. One exception to this pattern occurs in the vicinity of 44.58 N, 85.08 W, where a cluster of four inland wells exhibit relatively low d 2 H values compared to their neighbors.
Groundwater d values also exhibit significant spatial autocorrelation following removal of the observed trend with longitude (Moran's I ¼ 0.16, P ¼ 0.0001). We generated an interpolated prediction map for d using the IsoMAP geostatistical code and site longitude as a predictor variable (Fig. 4B). Cross-validation MAE and RMSE values for this model were 1.8% and 2.3%, respectively, representing a slightly larger fraction of the observed range of values than for the d 2 H model (15% and 19%). The interpolated d prediction map smoothes the variability observed at the individual monitoring sites, giving minimum and maximum values of 11.1% and 16.7% in the northeastern and southwestern quadrants of the Peninsula, respectively. Values are highest in the southwestern and northwestern parts of the study domain, and the map suggests a gradual decrease in groundwater d values toward the east and center of the Peninsula.
Because of the lower signal-to-noise ratio in the mapped d values we compare our prediction map with two spatially-explicit estimates of uncertainty in the map predictions (Fig. 5). The kriging interpolator gives an analytical expression for prediction uncertainty based on the spatial correlation structure of the observed data (Cressie 1993), and a map of the standard deviation of the kriging predictions shows a narrow range of values from 2.3% to 2.5%, with uncertainty decreasing with proximity to observation sites (Fig. 5A). Although these values represent a relatively large fraction of the total range of mapped values, there is no strong suggestion that specific features in the d map (Fig. 4B) are associated with regions of particularly poor predictions. Moreover, the values mapped in Fig. 5A provide an estimate of the uncertainty associated with a prediction at an individual point, but because the predictions are derived from a spatially continuous function, errors at different grid points are not indepen- v www.esajournals.org dent and the analytical uncertainties do not provide a robust basis for evaluating the strength of the predicted spatial patterns.
To assess the stability of the spatial patterns we conducted a full cross-validation of the interpolation process, creating 97 prediction surfaces based on independent model calibrations and interpolations using all possible subsets of 96 data. A map of the range (maximum À minimum) of these predictions at each grid point shows that at most locations predictions vary by less than 0.5% across all iterations, but also highlights several hot spots where predictions are more sensitive to individual observations (Fig. 5B). The highest range values occur in the vicinity of local 'outlier' observations. The presence of a small number (;6) of these outlier values in the data set is not surprising given the potential for secondary factors (e.g., evaporation of precipitation and infiltrating water, variation in recharge age), to produce groundwater d variation at local scales. Overall, the range values are higher in regions where mapped d values are highest and lowest, but the relatively small magnitude of the cross-validation ranges compared to the range of spatial variation in the d prediction map confirms that the mapped patterns are not particularly sensitive to individual observations.
Accepting that climatological variation in recycled water contributions to precipitation is the primary control on the large-scale patterns of d, the map values can be translated directly to values of the fraction of recycled water in precipitation ( f r ) using Eq. 2. The f r values are a linear transformation of the mapped d values and exhibit the same spatial pattern, with the highest fraction of recycled water in groundwater in the northwest and southwest corners of the Lower Peninsula and lower values characterizing the eastern and central parts of the study region. The f r values for the gridded groundwater map range from 3% to 12% using h ¼ 0.88 and from 5% to 21% using h ¼ 0.61 (see Eq. 1). In order to assess these values in the context of regional hydrological fluxes, we adopt the assumption that the groundwater represents an unbiased integration of local, annual precipitation at each grid cell and multiply the mapped f r values by annual precipitation amounts to obtain gridded estimates of the annual accumulation of recycled water in precipitation (Fig. 6). These values range from 22 to 109 mm/year for the h ¼ 0.88 case and from 38 to 198 mm/year assuming h ¼ 0.61, and exhibit a pattern of spatial variation similar to that of the d values. Summing the recycled water amounts across the study domain, we estimate that a total of between 4.8 and 8.7 km 3 of precipitation derived from recycled water falls within our study region annually. Of this, between 3.6 and 6.5 km 3 /yr falls within the Lake Michigan drainage basin and potentially contributes to stream and groundwater discharge returning to the Lake. For comparison, the annual recycled water estimates for the Lake Michigan basin equal approximately 0.8-0.14% of the volume of Lake Michigan and 9.8-17.8% of the total annual evaporative flux from the Lake surface (Grannemann et al. 2000).

DISCUSSION AND CONCLUSIONS
Isotope ratios of shallow groundwater in Michigan's Lower Peninsula vary over both local and regional spatial scales. Coherent, large-scale variation exists for d 2 H and d values, suggesting that these values reflect regional variation in the isotopic composition of local precipitation-derived recharge. In the case of d 2 H, the patterns largely parallel expected meridionally oriented gradients controlled by temperature patterns and the progressive distillation of water from northward-moving summer air masses (Rozanski et al. 1993). Our data set suggests that groundwater d 2 H values are several % lower within ;50 km of the Lake Michigan shoreline than they are further inland. Several factors could contribute to this pattern, including higher 2 H-depleted winterseason precipitation amounts in the vicinity of the Lake (Scott and Huff 1996) or enhanced contributions of relatively 2 H-depleted recycled Lake water to atmospheric vapor in this region. The data and map presented here provide a relatively high-resolution characterization of shallow groundwater isotope ratios that may be of use in characterizing water source d 2 H and d 18 O values in ecohydrological, ecological, paleoecological and forensic studies (Ehleringer et al. 2008, Wassenaar et al. 2009, Henne and Hu 2010, West et al. 2010, Kennedy et al. 2011. In the case of deuterium excess, the observed variation likely reflects regional differences in the amount of recycled water in precipitation, which can be quantified based on the model of Gat et al. (1994). Patterns of d variation suggest that the amount of recycled water in groundwater, and likely in precipitation, is highest close to the Lake Michigan coast, particularly in the northern and southern parts of the Peninsula. This pattern is consistent, in general, with well-documented patterns of lake effects on climate in this region (e.g., Scott and Huff 1996). In particular, the fall and winter season 'lake effect', which promotes precipitation from airmasses traversing Lake Michigan, is particularly strong along the west coast of the Peninsula, and may drive enhanced rainout of recycled water in near-coast regions (Scott and Huff 1996). As the dominant westerly air masses move further inland, lake effect precipitation decreases and moisture is added from transpiration and air mass mixing, diluting the recycled lake water contribution to vapor.
Localized recycled-water maxima in the northwestern and southwestern Peninsula are colocated with zones of particularly high fall and winter lake effect precipitation (Braham andDungey 1984, Scott andHuff 1996). In the northwest part of the Peninsula, this precipitation is dominantly snowfall, and may be driven by the relatively high topography of this region as well as the fact that northwesterly airmasses reaching the region are influenced by both Lake Michigan and Lake Superior. The highest rainfall amounts on the Peninsula occur within the southwestern zone of high recycled water values, and are related to high precipitation amounts in the late fall and early winter when temperature contrasts between Lake Michigan and overriding air masses are largest (Changnon and Jones 1972). In addition, the fetch of westerly and northwesterly air masses traversing the Lake is somewhat longer in this part of the Peninsula than at sites further north. The extent to which spatial variation in recycled water contributions to precipitation reflects regional variation in flux of water from the Lake to the atmosphere as v www.esajournals.org opposed to variation in the occurrence of 'lake effect' storms (i.e., the probability that the recycled vapor will be rained out) remains undetermined and will be important to evaluating projected future changes in regional hydroclimate. This issue could be addressed through spatially distributed monitoring of d in atmospheric vapor.
Our approach provides a unique opportunity to estimate the volumetric flux of recycled water in precipitation within a large part of the Lake Michigan drainage basin. These results suggest that between 10% and 18% of lake surface evaporation is currently re-precipitated within the leeward portion of the basin, although this value is sensitive to assumptions about the timing of aquifer recharge-if recharge occurs preferentially during storms or seasons with relatively high (low) recycled water fractions, the recycled water flux may be over (under) estimated. A significant fraction of this water is likely returned to the lake and represents a closed loop in the basinal water budget. Understanding the internal and external systems dynamics governing this loop, including the relative importance of over-lake evaporation rates, lake effect storms, and partitioning of land-surface water fluxes on recycled water return to the lake, will be critical to predicting the sensitivity of Lake Michigan water balance in response to climate and land use change. Such changes also have the potential to alter recycled-water 'subsidies' to terrestrial ecosystems downwind of the Lake, with potential long and short-term ecological and biogeochemical impacts. For example, previous work in Michigan demonstrates that lake-effect snow has strong impacts on ecological processes (e.g., soil development, nutrient cycling, soil water recharge) that affect ecosystem structure (Schaetzl 2002, Stottlemyer andToczydlowski 2006). Consequently, long-term ecological impacts of lake-effect precipitation are evident in the northern Lower Peninsula, where snowfall abundance has been linked to the spatial distribution of mesic forest types (Henne et al. 2007). The patterns of recycled-water subsidies demonstrated by our new data parallel these forest type gradients (e.g., as reflect by variation in abundance of key taxa including sugar maple, beech, and hemlock), and further suggest strong links between recycling and ecosystem structure within the study region.
Our work demonstrates that geospatial analysis of groundwater isotope data from a spatially distributed sampling network offers a tractable approach to the identification and quantification of regional land-atmosphere water fluxes. In this case, we were able to identify significant variation in the delivery of recycled water to terrestrial systems across the Lower Peninsula of Michigan, hypothesize relationships between these patterns and hydroclimatic mechanisms, and quantify the regionally integrated magnitude of this flux. As new measurement technologies (Worden et al. 2007, Gupta et al. 2009) and monitoring programs (Keller et al. 2008) increase the availability and density of spatially distributed atmospheric isotope data, similar analyses will provide a basis for the direct quantification of land-atmosphere fluxes without the uncertainties inherent to our extrapolation of atmospheric and land-surface processes from groundwater data. This work will provide a basis for the observation of regionally integrated recycled-water fluxes, providing new information on hydroclimatological dynamics and ecosystem water provisioning as well as previously unattainable constraints on models of the water cycle. SUPPLEMENT C++ code used to model and map water isotope ratios and accompanying documentation (Ecological Archives C003-005-S1).