Seasonal variation in tap water δ2H and δ18O isotopes reveals two tap water worlds

Stable isotope ratios of hydrogen and oxygen (δ2H and δ18O) in tap water provide important insights into the way that people interact with and manage the hydrological cycle. Understanding how these interactions vary through space and time allows for the management of these resources to be improved, and for isotope data to be useful in other disciplines. The seasonal variation of δ2H and δ18O in tap water within South Africa was assessed to identify municipalities that are supplied by seasonally invariant sources that have long residence periods, such as groundwater, and those supplied by sources that vary seasonally in a manner consistent with evapoconcentration, such as surface water—the proposed two tap water “worlds”. Doing so allows for the cost-effective spatial interpolation of δ2H and δ18O values that likely reflect that of groundwater, removing the residual error introduced by other sources that are dependent on discrete, isolated factors that cannot be spatially generalised. Applying the proposed disaggregation may also allow for the efficient identification of municipalities that are dependent on highly variable or depleted surface water resources, which are more likely to be vulnerable to climate and demographic changes.

Understanding the variation of stable hydrogen (H) and oxygen (O) isotopes in tap water has global application for several diverse fields 1 . In many contexts, the management of water resources is expected to become increasingly challenging given the predicted impacts of climate change, urbanization, economic development, and increasingly complex water resource supply systems and infrastructure 2,3 .
Municipal tap water is a primary interface between people and the hydrological system. Most municipalities are supplied by some combination of local surface and groundwater but may also rely on non-local sources through inter-basin water transfer schemes and processed water such as treated sewage effluent or desalinated seawater 4,5 . A theoretical understanding of the seasonal variability of stable isotopes in precipitation, surface water and groundwater ( Fig. 1) provides the foundation from which informed inferences can be made about the seasonal variability of stable isotopes in tap water. Some of these inferences have recently been corroborated 2,6-13 , but others remain untested.
Understanding the source of tap water supply is important for sustainable water resource management, as the associated risks, including those due to climate change, differ between sources 10 . Greater short-term variability in water input from precipitation, for example, would have a greater impact on hydrological systems that have a limited storage capacity, such as smaller surface water reservoirs, than those with longer residence times, such as groundwater aquifers or large surface water reservoirs. Furthermore, the monitoring of groundwater resources is critical to the sustainability of water supply, especially where groundwater resources are relied upon to buffer the effects of greater climate variability and where water demand has exceeded the supply capacity of traditional surface water resources. Information on the source of municipal tap water is not always available, reliable or independently verifiable, particularly in low-and middle-income countries. The available groundwater information in South Africa, which is one of the most data-rich countries in Southern Africa, for example, is inadequate for management purposes 14 and reliable records are only available for fewer than a third of the major surface water reservoirs 15 . Stable isotopes in tap water offer an approach to monitoring the dynamics and utilisation of critical water resources at a national scale.
The primary aim of this study is to test a cost-effective method of identifying and quantitatively characterizing the two tap water worlds-the first, groundwater with seasonally invariant δ 2 H and δ 18  www.nature.com/scientificreports/ tested at a national scale using tap water from South Africa as a case study. Successful disaggregation between groundwater and surface water sources would provide greater insight into the dynamics of the water supply system. It would enable the identification of municipalities that may be vulnerable to predicted climate and demographic changes, enable the prioritisation of adaptation interventions to increase the resilience of the water supply system, and it would provide a framework to directly and cost-effectively monitor the water resources on which municipalities are reliant.  8 .

Collection of samples.
Prior to analysis, each sample was filtered through a 0.2 µm syringe filter. Six injections of 1.9 µl were analysed for each sample, with the results of the first three injections discarded as conditioners. All samples were assessed for additional components that could affect the analysis following the methods of West et al. 8 . Isotopic ratios are expressed in ‰ as: . Surface water varies along the LEL due to evapoconcentration in the dry months and recharge in the wet months. Groundwater is comparatively invariant with isotopes in unconfined groundwater reflecting that of mean annual precipitation and confined groundwater differing from the local mean annual precipitation but rather reflecting that of the spatially or temporally distant recharge source (such as paleorecharge). The δ 2 H and δ 18  For drift correction and quality control, three in-house standards of known δ 2 H and δ 18 O values (1.6 and − 0.6; − 71.7 and − 10.2; and − 129.5 and − 17.27, respectively) where analysed after at most six consecutive samples were run. Tap water samples were calibrated to VSMOW using two in-house standards. The accuracy and precision of the calibration were assessed with a third quality control standard not used in the calibration. The mean (± one standard deviation) accuracy and precision of the isotope analyses across all batches were 0.4‰ ± 0.2‰ and 0.5‰ ± 0.2‰ for δ 2 H and 0.10‰ ± 0.05‰ and 0.09‰ ± 0.05‰ for δ 18 O, respectively. Spatial modelling. A Monte Carlo simulation permutation test for Moran's I statistic was calculated using 5,000 random permutations for the observed δ 2 H, δ 18 O and d-excess values in R version 3.6.3 (https ://cran.rproje ct.org/) 17 using the dismo 18 and spdep 19 packages. Universal kriging models were applied to interpolate the observed δ 2 H, δ 18 O and d-excess values and construct prediction isoscapes for each. The additional explanatory variables that were considered include the mean annual precipitation (MAP) and temperature (MAT), elevation (elev) and A-pan equivalent potential evaporation (PE), obtained from Schulze et al. 20 . The shortest Euclidean distance to the coast (tocoast) was calculated in R 17 using the rgdal 21 and sf 22 packages. Lastly, the predicted mean annual δ 2 H (MAH), δ 18 O (MAO) and d-excess (MAD) in precipitation [23][24][25][26] were also considered.
A linear model was constructed for each isotope ratio from the November and May tap water sampling campaigns using all five climatological and topographical variables (i.e., MAP, MAT, elev, apan and tocoast), the interaction between elev and tocoast (elev*tocoast) and between MAP and apan (MAP*apan), and the corresponding isotope ratio in precipitation (i.e., where δ 2 H in tap water is being predicted, MAH is considered www.nature.com/scientificreports/ as an explanatory variable but not MAO or MAD). The explanatory variables for the final model were selected by first assessing the stepwise effect on the AIC (Akaike information criterion) value and then the statistical significance of the variable coefficients. The kriging models were run in R 17 using the automap 27 , gstat 28 , sp 29 , and raster 30 packages, the and plotted using the tmap 31 package. Where the model uncertainty was greater than one standard deviation of the underlying training data, hatching was overlaid on the isoscapes to represent a lower confidence in the model predictions.
Quantification and interrogation of seasonal variability. Spatial modelling includes a degree of uncertainty and residual error in its predictions. To investigate the seasonal isotope variation mechanisms, therefore, further www.nature.com/scientificreports/ interrogation of the data was run on the observed isotope ratios for the 158 locations that provided samples for both campaigns rather than on the interpolated isoscapes that were predicted based on universal kriging. The seasonal variability in δ 2 H and δ 18 O of repeated sampling locations was classified based on the magnitude and gradient of variation (Δ) between May and November (Fig. 2) conservatively based on the characterization of surface water δ 2 H-δ 18 O evaporation slopes by Gibson et al. 32 . The thresholds applicable for a given locality differ based on numerous geographical, topographical and climatological factors 33 . The criteria applied in this study are therefore not globally generalisable and warrant further investigation. To allow for any variability resulting from sampling and observational errors 34,35 , a greater margin of variation was allowed for in the classification criteria than would be expected based solely on theoretical principles 33 . The classifications included seasonally invariant locations (Category A; |Δδ 2 H| < 5 and |Δδ 18 O| < 0.6), those that varied in a manner consistent with evaporation (Category B; 4 < �δ 2 H �δ 18 O <7), those that had an atypically shallow gradient of variation (Category C; . Each category was further distinguished by the direction of variation between May and November. Sample locations with isotope ratios in November that are more positive than in May are categorized as having a positive direction of variation (i.e., B p , C p and D p ) and vice versa for sample locations categorized as having a negative direction of variation (i.e., B n , C n and D n ).
The regional drift in the modelled isotope ratios is best explained by climatic and topographical variables, depending on the isotope ratio and the sampling period ( Table 1). The difference in the AIC values of the full and final models was relatively small but, to prevent over-fitting, variables that increased the model AIC were still removed. The highest R 2 value of the final models was 0.33, with p < 0.001, so a statistically significant amount of the null deviance was described in each case by the explanatory variables. The results of the Moran's I tests indicate that there is a significantly positive spatial autocorrelation in the isotope ratios (Table 2).

Quantifying seasonal variation. Seasonally invariant (Category A) samples made up 43% of all repeated
sample locations, with categories B, C and D making up 40%, 8% and 9% respectively (Fig. 5). The distributions of category A isotope ratios were consistent with the expectation for meteoric water (Fig. 6). Category B May samples were similar to those of category A, but reflected a greater degree of evapoconcentration and precipitation of recycled continental moisture in November than in May for the positive and negative subcategories, respectively (Fig. 6). Precipitation of recycled continental moisture was evidenced primarily by the more positive d-excess values, some of which were greater than 15‰ (Fig. 5). Evidence of recycled continental moisture precipitation was not observed in either category C or D and considerably negative d-excess values were only observed in some November C p samples (Fig. 6). Despite notable differences between May and November distri-  (Fig. 6). The variation in category D isotope ratios was not consistent with differences between the seasons in the degree of evapoconcentration or proportion of precipitation from recycled continental moisture. The distribution of the four categories of seasonal isotope variation varied spatially across the region (Fig. 7). Sample locations in the Northern Cape and North West provinces were dominated by category A variation, sample locations in the Western Cape, Eastern Cape, Free State and KwaZulu-Natal provinces were dominated by category B variation, Gauteng province consisted almost entirely of category D variation and the sample interpolation of variation categories. Isoscapes for the May (Fig. 8) and November (Fig. 9) δ 2 H, δ 18 O and d-excess values for category A and B sample locations were produced by universal kriging. The density and distribution of category C and D sample locations are too limited to allow for representative kriging of national scale isoscapes. The spatial variability of category A samples broadly reflected that of the predicted precipitation isotopes (Fig. S3). The δ 2 H and δ 18 O values for category A were more negative in the interior than on the coast; the eastern coast in KwaZulu-Natal province, in particular, had more positive δ 2 H values than the rest of the study area.
Several category B sample locations were discrete outliers, demonstrating various degrees of evapoconcentration. This was evident by 'bulls-eye' patterns of atypically positive δ 2 H and δ 18 O values and atypically negative d-excess values (Figs. 8 and 9). These discrete evapoconcentrated outliers were more obvious in the November than the May sampling campaign.
The differences in the spatial coherence of samples from the two categories were evident in the residual error of the universal kriging models (Figs. S7 and S8). The coherence of category A samples were reflected in the low residual error, even in regions with a low density of observations such as the Eastern Cape province (see Fig. 7 for the distribution of category A sample locations). Conversely, the lack of coherence of category B samples was reflected by the high residual error at short distances away from the observed sampling locations. Model uncertainty greater than one standard deviation of the underlying data was especially evident in the November category B isoscapes (Fig. 8). This indicates that there is less confidence in the interpolation of observed isotope ratios beyond the sampling location for this category. www.nature.com/scientificreports/

Discussion
We demonstrate that remote, citizen-sampled seasonal variability of tap water isotopes is potentially an efficient and cost-effective method of identifying and quantitatively characterizing groundwater and surface water dependent tap water sources-the two tap water "worlds". Seasonal variability in stable water isotope ratios, as defined by the magnitude, gradient and direction of change between May and November, provides a powerful diagnostic tool to help identify the primary source of tap water supply. Sampling locations that exhibited a low magnitude of seasonal variability (Category A; |Δδ 2 H| < 5‰ and |Δδ 18 O| < 0.6‰) are consistent with the expectations for groundwater 8 . Sampling locations that exhibited seasonal variability that approximate local evaporation (Category B; 4 < �δ 2 H �δ 18 O < 7) are consistent with the expectations for surface water 36 . After disaggregating between surface and groundwater reliant municipalities, considerable distinctions are noted in the spatial variability of δ 2 H and δ 18 O (Figs. 8 and 9). This has direct implications on the application and utility of stable isotopes in tap water for water resource management. Surface-coupled water sources are exposed to risks that are distinct from that of groundwater reliant sources; for example, they are more impacted by silting up due to upstream erosion processes. They are also at greater risk of increased evaporative losses due to climate change. Groundwater sources, by contrast, are more at risk of saltwater intrusion and from dissolving www.nature.com/scientificreports/ salt from the substrate. Furthermore, this has implications for how stable isotopes are applied in fields such as human (ground and surface water use) and wildlife (surface only) forensics, which rely on the predictability of spatiotemporal variation in stable isotope ratios 37 . The climatological and topographical explanatory variables used in the kriging models accounted for less variation in the isotope ratios from the November campaign than those of the May campaign (Table 1). This may be a result of the greater degree of evapoconcentration that was observed, on average, in November than in May samples (Fig. 3). The degree of evapoconcentration is dependent, in part, on the surface area to volume ratio of the evaporating water body with greater ratios resulting in a greater proportion of evaporate to residual water than reservoirs with smaller ratios. The spatial distribution of evapoconcentration is therefore not primarily dependent on climatological or topographical variables that vary continuously, but on discrete reservoirspecific characteristics. Therefore, the more that the spatial variation in δ 2 H and δ 18 O values are determined by evapoconcentration, the less the predictive power of climatological and topographical variables or spatial autocorrelation will be. This is also reflected in the greater standard deviations of the category B kriging models for November compared to May in areas with a lower sampling resolution (Figs. S7 and S8). The extrapolation of δ 2 H and δ 18 O values beyond the location of sampling is unreliable for this category (Fig. 8) as the mechanism www.nature.com/scientificreports/ of variation is dependent on discrete factors such as the properties of the reservoir source. Linear interpolation based on environmental factors and spatial autocorrelation is therefore inappropriate. The sample locations that demonstrated medium to high seasonal variability (Δδ 2 H| > 5‰ and |Δδ 18 O| > 0.6‰) have several possible interpretations. The source of tap water supply could entirely or partly include surface water reservoirs that experience a greater degree of evapoconcentration in the dry than the wet season, or the source could entirely or partly include non-local water that has stable isotope ratios distinct from that of local water. Non-local water could be distant in space (i.e., inter-basin water transfer from a different catchment area or desalinated ocean water) or distant in time (i.e., fossil groundwater that was recharged under different climatic conditions). The gradient of change in water stable isotopes �δ 2 H �δ 18 O can help distinguish between the possible scenarios.
The seasonal variation from local evapoconcentration can only result in isotope ratio changes along the local evaporation line (LEL) in the positive direction 36 . If a given municipality blends some ratio of local ground and surface water, the gradient of change in the isotope ratios will still not deviate from the LEL if all local water sources are recharged by local precipitation. Constrained groundwater aquifers that were recharged under  www.nature.com/scientificreports/ different conditions may have isotope ratios that differ from the local mean, in which case blending or sourceswitching between the local and non-local water will result in deviations from the LEL. The same is true for inter-basin water transfer (IBWT) schemes 12 . Figure 10 presents some hypothetical seasonal variation scenarios following the blending of isotopically distinct local and non-local water. The assumption is that there is a greater degree of evapoconcentration in the local water reservoirs that require augmentation given that they would be more depleted. If the non-local water source is a constrained groundwater aquifer, rather than from IBWT, then similar patterns in the seasonal variation would occur-the only difference being that the non-local source would experience no evapoconcentration and would be isotopically invariant. The direction of change depends on whether winter or summer is the dry season for the given sample location. Sites in summer rainfall regions would be expected to vary in the opposite direction to those in winter rainfall regions. Sampling locations in the Gauteng province offer an empirical example of an IBWT scheme 4 .
The water requirements of the urban municipalities and industrial sectors in Gauteng are met by a large IBWT scheme sourced from Katse Dam in Lesotho 4 . Two expectations are provided, based on the principles outlined above. The climatic and topographical dynamics between the Gauteng and Lesotho catchments suggest that the scenarios illustrated in Fig. 10a best reflect that of this case study. The degree to which the Gauteng municipalities are dependent on water supply from Lesotho in the dry season suggests that the scenario illustrated with the green line (positive gradient) better reflects this case study than that of the purple line (negative gradient) in Fig. 10a. These two expectations are supported by the direction and the gradient of variation, respectively, in tap water δ 2 H and δ 18 O values between May and November (Fig. 11).
Future studies can improve the interpretation of seasonally variable tap water isotopes by quantifying the spatial and seasonal variability in the LEL gradient. Additionally, further comparisons between the nationalscale groundwater isotope ratio predictions based on tap water samples and catchment-scale groundwater isotope data would be valuable to determine the accuracy at scales relevant to local water resource management. Comparisons between the predicted groundwater isotope ratios and that of local modern precipitation would also contribute to determining whether or not the groundwater being extracted was likely recharged recently or is non-renewable 12 . The use of other geochemical tracers may further distinguish between modern and fossil groundwater aquifers 38,39 .
The seasonal variability of δ 2 H, δ 18 O and d-excess values in tap water aid in to their application to urban hydrology. In particular, it has allowed for the identification of the two tap water "worlds"-municipalities that likely rely on groundwater and surface water dependent tap water sources. This has several implications for fields such as water resource management, ecology and forensics 1 . This study has demonstrated that stable H and O isotopes in tap water have the potential to complement traditional water resource management practices. The degree of observed temporal variability suggests that future studies of stable isotopes in tap water need, at a minimum, to incorporate seasonal sampling in order to understand the variable source and transport dynamics of a given location.