Climatic and Topographic Control of the Stable Isotope Values of Rivers on the South Island of New Zealand

We show that climate and topography control the spatial distribution of stable isotope values on the South Island of New Zealand, based on a spatially dense (n = 193) river isotopic survey. Our data show a δ18O minimum in isotope values east of the Southern Alps that demonstrates topographically driven continentality associated with the Southern Alps, which intersect the prevailing, moisture‐laden westerlies. Our data define a South Island surface water line of δ2H = 8.17 (±0.26) × δ18O + 10.57 (±2.04), which is identical within 95% confidence intervals to the global and New Zealand meteoric water lines established from monthly precipitation samples. The observed river δ18O values are strongly correlated with annual temperature range and winter temperature. Strongest correlations are between δ18O and mean minimum winter temperatures (r > 0.7 for June, July, August), with gradients of 0.58–0.66‰ /°C. Based on a multiple regression analysis of δ18O against climate data, we present a river δ18O model and isoscape that demonstrate the control of continentality and moisture source on New Zealand surface water isotope spatial patterns. Model validation against previously published river samples shows skill in predicting river δ18O values (root‐mean‐square error = 0.83), confirming that the spatial variations in river δ18O (and δ2H) are robust to sampling period and reflect continental, precipitation source and temperature effects. Our data suggest that oxygen or hydrogen isotope paleoclimate proxies derived from rivers or open‐system lakes on the South Island should be sensitive to winter temperature.


Stable Isotopes of Environmental Waters
Temporal and spatial variations in precipitation isotope values encode climatic processes that reveal atmospheric dynamics, influence of orographic barriers, and temperature controls. As such, δ 18 O and δ 2 H values of meteoric water are strong candidates to be considered as essential climatic variables (Bojinski et al., 2014) that have high utility for investigating modern and past climate. In mid-to high-latitude regions, the δ 18 O and δ 2 H values of meteoric waters are strongly related to temperature, which controls both the equilibrium fractionation between vapor and condensate (Fricke & O'Neill, 1999;Rozanski et al., 1993), and the amount of heavy isotope distillation associated with air mass rainout along an ocean to continent cooling path. Stable isotope values vary with air mass evolution within storms and with orographic rainout (Lachniet et al., 2016;Poage & Chamberlain, 2001;Rowley & Garzione, 2007;Winnick et al., 2014). Furthermore, stable isotopic changes relating to moisture sources arriving at a location are strongly related to subsequent weather system history and regional geography. The temperature control on meteoric water isotope Abstract We show that climate and topography control the spatial distribution of stable isotope values on the South Island of New Zealand, based on a spatially dense (n = 193) river isotopic survey. Our data show a δ 18 O minimum in isotope values east of the Southern Alps that demonstrates topographically driven continentality associated with the Southern Alps, which intersect the prevailing, moisture-laden westerlies. Our data define a South Island surface water line of δ 2 H = 8.17 (±0.26) × δ 18 O + 10.57 (±2.04), which is identical within 95% confidence intervals to the global and New Zealand meteoric water lines established from monthly precipitation samples. The observed river δ 18 O values are strongly correlated with annual temperature range and winter temperature. Strongest correlations are between δ 18 O and mean minimum winter temperatures (r > 0.7 for June, July, August), with gradients of 0.58-0.66‰ /°C. Based on a multiple regression analysis of δ 18 O against climate data, we present a river δ 18 O model and isoscape that demonstrate the control of continentality and moisture source on New Zealand surface water isotope spatial patterns. Model validation against previously published river samples shows skill in predicting river δ 18 O values (root-mean-square error = 0.83), confirming that the spatial variations in river δ 18 O (and δ 2 H) are robust to sampling period and reflect continental, precipitation source and temperature effects. Our data suggest that oxygen or hydrogen isotope paleoclimate proxies derived from rivers or open-system lakes on the South Island should be sensitive to winter temperature.

Plain Language Summary
We investigated the spatial variations in oxygen and hydrogen stable isotopes in rivers on the South Island of New Zealand and show that they are strongly controlled by the presence of the Southern Alps. Indicators of continental climates, such as seasonal temperature range and temperature of the coldest months are strongly related to measured oxygen isotope values of river waters. We show that winter temperature and continentality are the dominant control on New Zealand South Island river oxygen isotope values, and they result in minimum values in the lee of the highest portion of the Southern Alps. We also show that river deuterium excess values are controlled by Pacific Ocean versus Tasman Sea moisture sources.
values has led to their wide use in geologic materials, such as lake sediments (Leng & Marshall, 2004), ice cores (Brook & Buizert, 2018), soil carbonates (Oerter & Amundson, 2016), and cave calcite deposits (Lachniet, 2009;Nava-Fernandez et al., 2020;Whittaker et al., 2011;Wong & Breecker, 2015) as temperature-sensitive paleoclimate proxies. Stable isotopes are also routinely measured in biologic archives, for example in tree ring cellulose (Lorrey et al., 2016;McCarroll & Loader, 2004), to understand modern processes governing plant water use. Further, understanding the controls on modern environmental water δ 18 O and δ 2 H can be used to better interpret past climatic processes. As such, constraining the quantitative relationships between environmental water δ 18 O and δ 2 H values relative to a wide range of physical processes is a key goal of modern isotope climatology. Spatial analysis of meteoric water isotope values also allows the production of interpolated maps of the δ 18 O and δ 2 H composition of waters called isoscapes (Bowen & Wilkinson, 2002;West et al., 2010). A meteoric precipitation isoscape has been completed for New Zealand based on a network of 51 collection stations on the North and South Islands (Baisden et al., 2016). These data can be used to constrain moisture sources and physiographic and climatic controls on precipitation (Baisden et al., 2016;Frew et al., 2011;Rogers et al., 2012). The monthly precipitation isotope data were explained by altitude, temperature, and precipitation amount. However, continuous and long-term isotope precipitation collection is costly and time intensive (Blackstock, 2011), and few regions contain a sufficient density of precipitation isotopic monitoring stations to resolve spatial variability. This problem is particularly acute in regions of high topographic diversity like New Zealand. As a result, studies based only on meteoric precipitation monitoring may miss important geographic variability in isotope values. A geographic survey of environmental waters from other sources, such as rivers, can permit a more comprehensive investigation of spatial variations.
The spatial variations of stable isotopes have powerful geological applications. δ 18 O contained in natural geological archives have been used to constrain mountain uplift and the development of negative isotopic anomalies in the lee of high mountain ranges (Mulch, 2016;Poage & Chamberlain, 2001;Rowley & Garzione, 2007). In New Zealand, for example, a 5-6‰ δ 18 O decrease in authigenic kaolinites in the lee of the high Southern Alps of the South Island was used to infer mountain uplift and the establishment of a continental rain shadow in the early Pliocene (Chamberlain & Poage, 2000;Chamberlain et al., 1999). The modern calibration data set of water δ 18 O values used to infer isotopic gradients for that study was small (n = 24), lacked geographic coordinates for independent verification, and included lakes which may have been evaporatively enriched in heavy isotopes and therefore potentially mischaracterized the magnitude of the isotopic gradients. A more comprehensive suite of river water isotope values was presented by Stewart et al. (1983), who measured δ 2 H values in 61 rivers on the South Island, with most samples located near the coasts. Additionally, Kerr et al. (2015) sampled 29 closely spaced sites in a transect from the Tasman Sea to near the alpine divide. These three studies indicate that river isotope values vary significantly across the South Island, but large areas remain unsampled. To improve the utility of spatial δ 18 O variations for climate and environmental research applications, more detailed studies based on larger data sets are useful to evaluate the climatic control on environmental water isotope values.
Here, we combine 193 measurements of H and O stable isotopes from rivers on the South Island of New Zealand collected mostly in November of 2016 to constrain the climatic controls on river water δ 18 O values (Table S1). This data set is more than three times larger than previously published data sets and spans a greater topographic range and representation in South Island climate districts. We correlated the river isotope data with climate data sets to investigate the climatic and physiographic controls on river δ 18 O and δ 2 H values of our samples. We then built a multiple regression model for the South Island that serves as a baseline from which to compare future river sampling campaigns. The objective of this study is to better constrain the climatic controls on the spatial variation in river water δ 18 O, δ 2 H, and deuterium excess values. In addition, we develop a spatial representation of river isotope variations that can be utilized in future applications for isotope climatology, geo-and bio-forensics, and water resources management, and to guide site selection for future paleoclimatic and paleotopographic analysis. Our conclusions have direct bearing on the interpretation of isotopic paleoclimate records, and thus serve as a modern constraint on how past climate variations could be encoded in variations of oxygen and deuterium isotopes. This paper presents an open-access surface water database available to the scientific community for use in more detailed investigations at finer spatial and temporal scales.

River Water Stable Isotopes
River waters inherently capture spatio-temporally averaged processes related to the complex interplay of geology, geography, climate, and weather. This means river isotopes can be influenced by mean residence time of water in a catchment, changing and variable surface-groundwater proportionalities, contrasting inter-and intra-season meteoric water contributions, snowmelt, and reach-specific baseflow changes. River water isotopic compositions are related to the δ 18 O and δ 2 H values of meteoric precipitation that falls within a catchment (Kerr et al., 2015;Lachniet et al., 2016;Mulch, 2016;Stewart et al., 1983). As a result, rivers commonly experience seasonal variations in isotopic values that respond to the changes in isotopic values of meteoric precipitation that falls into their catchments but are dampened and delayed in time. Rivers are thus biased to represent the stable isotopic values of the meteoric precipitation and groundwater entering the catchments prior to the period of sampling, and these values may change over time. Ideally, a long-term sampling campaign on each river can establish mean δ 18 O values and estimates of isotopic seasonality over time, but such a study becomes difficult to apply for the 193 locations at which we sampled. Fortunately, previous studies (Stewart et al., 1983) have shown that the seasonal variations in New Zealand stream water isotopes are small. For example, δ 2 H variations of up to 8-10‰ were observed for catchments with significant snowmelt contribution; when converted to δ 18 O using the global meteoric water line (GMWL), this represents seasonal variations on the order of 1.0-1.2‰. However, other rivers studied by Stewart et al. (1983) have seasonal δ 2 H variations of only 1-5‰, corresponding to δ 18 O variations of ∼0.2-0.6‰.
While much attention has been given to the temporal variations in precipitation stable isotope values, the spatial variation in river stable isotope values may also reveal important controls from climatic and geographic processes. Geographic ranges in δ 18 O values are commonly large relative to the amount of seasonal variation, such that spatial patterns of δ 18 O are robust to the time of sampling. For example, a river water δ 18 O isoscape for Alaska showed spatial changes of up to 14‰, characterized by a decrease in river δ 18 O values with increasing distance from the Gulf of Alaska (Lachniet et al., 2016), with associated minima over mountain ranges associated with the lower condensation temperatures over high terrain. A multiple regression model for that region against climatic and physiographic data explained river δ 18 O values to better than ±2.0‰ and was consistent with winter precipitation stable isotope values at most sites where meteoric waters were sampled. Similar strong correlations between river δ 18 O values and climatic data was evident in the spatial variations in Central American rivers (Lachniet & Patterson, 2006, 2009. A decrease in δ 18 O values of rivers with increasing distance from the Caribbean Sea was interpreted to reflect the progressive distillation of heavy isotopes as air masses rained out upon traversing Central America (Lachniet & Patterson, 2009). A statistical model of river δ 18 O values based on climate and physiographic data for data from Guatemala and Belize explained river δ 18 O values to better than ±1.0‰. The coherent regional-scale river δ 18 O variations in both Alaska and Central America demonstrate that the noise associated with local watershed hydrologic processes are of lesser strength than the signal embedded in the spatial variations. The strong correlations between spatial δ 18 O values and climate variables in these studies demonstrated the utility of river water isotope values to encode climatic information, and validation of the statistical models with independent data allowed prediction within the ±1-2‰ range for river δ 18 O, which is considerably narrower than the range in spatial δ 18 O variation. This point should be kept in mind when interpreting the data in this study.
Seasonally resolved single site and regional snapshot surface water sampling campaigns provide complementary data to evaluate climate-isotope linkages. The climatic controls on water δ 18 O and δ 2 H values are also required to confidently interpret groundwater flow dynamics (Taylor et al., 1989), and they underpin paleoclimatic and/or paleoenvironmental proxy record interpretations from caves (Hellstrom et al., 1998;Lorrey et al., 2008Lorrey et al., , 2020Nava-Fernandez et al., 2020;Whittaker et al., 2011;Williams, 1996) and other isotopic registers like carbonates, molecular biomarkers, paleosols, and biogenic apatite (Browne et al., 2017;Hinojosa et al., 2019;Lorrey & Bostock, 2017;Mulch, 2016). As such, analysis of closely spaced networks of river water isotopes improve our understanding of the geographic variations in stable isotopes and are well suited for answering questions about the climatic and physiographic forcing on spatial isotopic variations.
River isotope data should be most representative of spatial variability when local seasonal variations are small. Our technique assumes that the relative spatial variations in river water δ 18 O values are related to primary climate processes, even if the absolute values may change over time because of seasonality or climate change. The available data for New Zealand (see above) indicate limited seasonal isotopic variability (Stewart et al., 1983), particularly for low elevation catchments without significant snowmelt contributions. The potential for variable meteoric source water contributions and water ages to a catchment can also be minimized by sampling small drainage basins that contain river water with precipitation sources that are located close to the sample site. Furthermore, rivers within small catchments better represent local climatic conditions than larger rivers with distant headwaters and sources (Kerr et al., 2015;Stewart et al., 1983). Spatial analysis of river water isotopes also supports sampling strategies that account for isotopic idiosyncrasies arising from significant topographic complexity, which occurs across the varied terrain of the South Island. Inevitably, each river-specific isotopic survey has the potential to add to the historic accountancy of surface water resources, which can be used as baselines from which we can establish spatial variability that arises naturally, compare on-going observations, and evaluate scenarios of future climate changes.

Climate of New Zealand's South Island
New Zealand's South Island has a midlatitude temperate and mostly maritime climate. It is positioned astride the prevailing southern hemisphere westerlies which interact with high topographic relief of >3,000 m to produce abundant precipitation on western slopes and upland terrain (Macara, 2015;Salinger, 1980). The Southern Alps, which began forming ∼5 Ma during the Miocene (Batt et al., 2004;Lorrey & Bostock, 2017), present a nearly continuous topographic barrier to incident atmospheric flow that forces air masses to rise, cool, and precipitate high amounts of rain and snow upon high terrain. Westerly air masses traversing the central Southern Alps exhibit only minimal flow deflection because of the uninterrupted high barrier of the Southern Alps (Wheeler & Galewsky, 2017). As a result, the interior South Island east of the main divide is characterized by a semi-arid continental climate that is significantly drier than the west coast and main divide (Macara, 2015). For example, the mean annual precipitation is <1 m in many areas east of the main divide compared to 12 m or more along the alpine slope facing the Tasman Sea and in the high Southern Alps (Henderson & Thompson, 1999). Anomalously high precipitation amounts on the South Island occur when there is a strong zonal pressure gradient (high in the north, low in the south), and stronger westerly winds which guide cyclonic disturbances to the east across the South Island (Salinger, 1980;Ummenhofer & England, 2007). The east coast receives some moisture from easterly sources (Macara, 2015;Stewart et al., 1983), and southwesterly winds and associated troughs deliver abundant precipitation to the south coast of the Southern Island (Macara, 2013).

Methods
We collected 193 water samples, mostly from rivers at sites across the South Island of New Zealand (Figure 1, Table S1), with most (n = 183) taken on a 5,550-km driving campaign November 13-29, 2016. Samples were collected from a fast-flowing portion of the stream using a purpose-built 3.7 l polyvinyl chloride cylinder with 30 m of polypropylene attached to a bail (aka bucket and rope), sterilized after each sampling with a concentrated bleach solution, and transferred to 50 ml HDPE plastic bottles with zero air headspace and the caps taped in place to prevent exchange with outside atmosphere. Samples ranged from 46.97°S to 40.73°S latitude and 167.1°E to 174.1°E longitude at altitudes from 1 to 826 m above mean sea level. Geographic coordinates for each sampling location (commonly taken on bridges spanning the rivers) were determined using a hand-held GPS relative to the WGS84 datum in decimal degrees. Sample locations included every physiographic region of the South Island, including all coastlines with transitable paved roads, with the exception of the temporarily inaccessible stretch between the locations of Kaikoura and Domett on the northeast coast resulting from road closures after the November 14, 2016 magnitude 7.8 Kaikoura earthquake (Cesca et al., 2017). We primarily targeted streams with small drainage areas that are likely to contain a local isotopic fingerprint close to the sample location. As a result, many large rivers on the South Island (e.g., the Clutha, Mataura, Waitaki rivers) were omitted from our collection. A few large rivers were analyzed but not included in the statistical analyses (e.g., the Waimakariri, Ashley, Hurunui, Awatere, and Wairau rivers). The subset of data used for statistical analyses was n = 186 samples. The average spacing between streams was approximately 34 km per sample, and our sample density was chosen to best balance the trade-off between sample density, spatial extent, and time dedicated to the sampling campaign. To assess the range of seasonal δ 18 O variations, we used data on 24 streams in the southern South Island, at five different dates spanning a year (Table S2), but not included in our main data set. Our data represent a snapshot of river δ 18 O values at the time of sampling and the focus of this study is to constrain the relative spatial variations. Although there are seasonal variations in river δ 18 O values, we expect, however, that the relative spatial variations in river isotopes are largely robust over time.
Water samples were analyzed for δ 18 O and δ 2 H using a ThermoElectron high temperature thermal conversion elemental analyzer (TC/EA) coupled to a Delta V Plus mass spectrometer with a ConFlo III interface in the Las Vegas Isotope Science Lab at the University of Nevada Las Vegas. Samples were reacted at 1400 °C in a glassy carbon tube packed with glassy carbon granules in a stream of ultra-high purity Helium and resulting CO and H 2 were separated on a GC column. Each sample from a 2 ml vial with zero headspace was injected into the TC/EA with a GC PAL autosampler six times, and δ 18 O and δ 2 H values estimated using CO and H 2 reference gases. The measured δ values were determined on the last three to five injections to avoid any memory effect, which tended to be small when geographically adjacent samples were analyzed sequentially. Sample analyses were completed in 2017 within 8 months of collection. An offline correction was applied using internal water standards that have values defined on a scale in which δ 18 O and δ 2 H of SMOW are both 0‰ and SLAP (Standard Light Antarctic Precipitation) are exactly −55.5 and −428‰ VSMOW, respectively (Thomas et al., 1996). Each run contained at least two internal standards with bracketing δ values (from either pair of LVOW1_2017 and DI-1_2017, or US Geologic Survey [USGS]-45 and USGS-47), and one internal standard (KIWI-2017) was processed and treated as an unknown as a check for internal reproducibility as defined by an overall standard deviation. Replicate measurements for the KIWI-2017 standard (n = 31) returned values of −4.57 ± 0.14 and −26.8 ± 1.9‰ for δ 18 O and δ 2 H, respectively. Average values and precisions for LVOW1_2017 were δ 18 O = 1.18 ± 0.13 and δ 2 H = −41.8 ± 2.38‰ (n = 18), and for DI-1_2017 were δ 18 O = −11.73 ± 0.17 and δ 2 H = −95.63 ± 1.53‰ (n = 18) VSMOW, respectively.
To investigate the potential climatic controls on the surface water isotope values, we correlated river δ 18 O to 112 climatic variables, many of them monthly, including from the WorldClim v2 and Bioclim databases at 30 arc-second spatial resolution (Fick & Hijmans, 2017), and the Global Aridity Index and potential evapo-Transpiration (PET) Climate Database v2 (Trabucco & Zomer, 2018), using least squares regression. Although our samples were collected mostly in November, we correlated the measured δ 18 O values to mean monthly climate parameters to test for the strongest potential controls. We also correlated our surface water stable isotope data against an additional five physiographic variables including sample latitude, longitude, distance from the coast (in decimal degrees), sample altitude, and catchment altitude (for a subset of rivers; n = 157) for a total of 117 different climatic and physiographic variables. The interpolated climate data, results of the correlation matrix, including Pearson r, p-values, slopes, and intercepts, are shown in Table S3.
The WorldClim database is available open-access to maximize utility to the research community. World-Clim version 2 represents monthly averaged climate data for 1970-2000 years CE from spatially interpolated thin plate splines and is improved in areas of large spatial gradients and rain shadows. For New Zealand, the WorldClim data is covered by 598 stations (Figure 1), and shows a high density of stations covering most of the South Island, with the exception of Fiordland (Fick & Hijmans, 2017). We also used the 1981-2010 NIWA precipitation coverages (https://niwa.co.nz/climate/research-projects/national-and-regional-climate-maps). These data were made using a trivariate thin plate smoothing spline that included latitude, longitude, and an expert-guided map of mean annual precipitation (Tait & Zheng, 2007), and thus are likely to be a more accurate representation of precipitation amounts and patterns for New Zealand. From the NIWA coverages, we used monthly average precipitation, mean annual precipitation, and seasonal precipitation totals (DJF, MAM, JJA, SON) and calculated precipitation of the wettest and driest months and quarters. The catchment parameters were derived from a USGS Earth Explorer digital elevation model (DEM) of New Zealand at 30-m spatial resolution, via mapping of drainage networks to delineate the upslope catchment of the water sample locations. Several built-in Spatial Analyst algorithms of ArcGIS 10.3.1 toolbox and a Python program was developed and executed to automate the delineation process. The DEM reconditioning was performed where elevation of any sink-hole pixels was raised to match lowest neighboring pixels.
The derived stream network provides a detailed surface flow pattern that can be used to delineate sub-catchments and upstream hypsometric mean, maximum, as well as catchment relief above the sample collection site. The stream network of the South Island was mapped in three steps. In the first step, a flow direction map was calculated from the reconditioned DEM based on eight-pixel neighborhood where the direction of flow at a pixel is toward the lowest neighbor. In the second step, a flow accumulation map was created and thresholded. In a flow accumulation map, a pixel value is the number of upslope contributing pixels. A threshold was applied to flow accumulation to identify locations where sheet flow becomes concentrated flow, that is, all pixels larger than threshold belong to streams. The threshold was chosen iteratively to ensure that the streams of water sample locations were delineated. Lastly, in the third step, a stream ordering algorithm was applied to group the pixels that form a single link of the stream network. Our approach resulted in 157 river catchments that could be matched with GPS coordinates; many of the other sampled streams were very small and were not captured because the resolution of the DEM was too low. In this research, the threshold was chosen through iterative adjustments to ensure streams for the maximum number of water sample locations. The water sample locations were snapped to the closest streams and their catchments were delineated. The catchment parameters for each sampled were mean, minimum, and maximum altitudes and upstream catchment area. To determine the altitude effect of stable isotopes, we used the mean catchment altitude instead of sample collection altitude.

Spatial Variations in River Water δ 18 O and δ 2 H Values
We observe clear spatial variations in stable isotope values of river waters. δ 18 O values were high on the Tasman Sea slope of the South Island (Figure 1), with a range between −4 and −7‰ for nonglacierized basins (e.g., Punakaiki and Nile Rivers). Glacierized catchments on the Tasman slope include an integration of ice and snow melt with rainfall up to the highest altitudes (>3,000 m) of the Southern Alps, and these rivers have lower δ 18 O values (e.g., Fox and Waiho rivers with δ 18 O values of −9.7‰ and −9.1‰, respectively) relative to nearby nonglacierized catchments. However, the lowest river δ 18 O values were found southeast of the Mount Aspiring and Aoraki/Mount Cook National Parks, where the continental divide is defined by peaks with altitudes higher than 2200 m. The stream with the lowest δ 18 O value we measured is the nonglacierized Luggate Creek (δ 18 O = −12.8‰), which has a δ 18 O value significantly lower than the glacially sourced Tasman River (δ 18 O = −9.6‰) with headwaters originating on Aoraki/Mount Cook at 3,724 m. In broad terms, the river δ 18 O values decrease from the Tasman Sea to east of the divide, then increase from just east of the Southern Alps main divide toward the Pacific coastline where δ 18 O is mostly between −7‰ and −10‰, with the exception of higher (∼−5‰ to −7‰) values around the southeast coast ( Figure 1). These observations indicate large spatial variations on the stable isotope values of South Island river waters. The highest South Island δ 18 O values occur along the upwind (western) slopes of catchments that drain toward the Tasman Sea, while the lowest values are from interior locations in the lee of the highest Southern Alps peaks. Intermediate to high δ 18 O values are also located along the Pacific slope.
Our data define a surface water line (Figure 2) for the South Island of δ 2 H = 8.17 (±0.26) × δ 18 O + 10.57 (±2.04), which is identical within 95% confidence interval uncertainties to the GMWL of δ 2 H = 8.17 × δ 18 O + 10.35 (Rozanski et al., 1993), and to the New Zealand meteoric water line δ 2 H = 7.92 (±0.05) × δ 18 O + 11.07 (±0.31) defined by (Frew et al., 2011). The δ 18 O/δ 2 H slope of 8.16 and intercept of 10.57 demonstrates that the river waters have not been significantly affected by evaporation, which, if a strong influence, would lower the slope and intercept away from GMWL values. The river waters show a smaller δ 18 O range than that observed in monthly precipitation values spanning 3 years (Frew et al., 2011), and our data range (  ure S2) compares better with July-August (mid-to-late winter) precipitation (Baisden et al., 2016). We attribute this to our sampling campaign at the end of the spring melt season (November 2016), when rivers contained snowmelt that fell as winter precipitation in the preceding months.

Climatic and Physiographic Controls on River δ 18 O Values
We observe strong correlations (Table S3)  The temperature gradients and 95% confidence intervals of δ 18 O versus T min range from 0.58 ± 0.07‰ per °C (July) to 0.90 ± 0.18‰ per °C (March). The gradients for T max are largest during late autumn to mid-winter (MJJ) with slopes of 0.42 ± 0.16 to 0.53 ± 0.13‰ per °C, and for T avg range from 0.72 ± 0.11 to 0.78 ± 0.14‰ per °C across mid-autumn through early spring (AMJJAS). Correlations to summer temperatures are weaker. Among the nontemperature variables, δ 18 O showed a strong inverse correlation to December potential evapotranspiration (r = −0.75), and PET and T max were the only ones that showed stronger significant correlations to summer climate data. Several variables had significant correlations to δ 18 O but mostly lacked a clear seasonal cycle: in order of decreasing strength of correlation, δ 18 O was correlated to vapor pressure (r ≥ 0.62), wind velocity (r ≥ 0.46), then precipitation (r = 0.32-0.40, NIWA data). Solar radiation had weak correlations to δ 18 O. In the BioClim data set, the strongest correlations were with temperature seasonality (Bio04, r = −0.75), minimum temperature of the coldest month (Bio06, r = 0.78), and temperature annual range (Bio07, r = −0.77) of which each parameter exhibits the most extreme values in the continental interior of the South Island. In contrast, correlations to precipitation parameters are lower than those for temperature, with the strongest correlations to precipitation of the driest month (r = 0.32, NIWA data) and precipitation of the wettest quarter (r = 0.37, NIWA data). The strongest correlation with the monthly precipitation data was for June precipitation (r = 0.40), which is relatively weak in comparison to the winter temperature variables.
Many of the climate variables are collinear and are related to the continental nature of the central South Island east of the Southern Alps. For example, distance from the coast is a strong predictor of minimum winter temperatures (e.g., r = −0.83 between July T min and distance from the coast), temperature seasonality (r = 0.85), and mean temperature of the coldest quarter (r = −0.74; see Table S3). In addition, distance from the coast is also a predictor of December potential evapotranspiration (r = 0.76), and several other climate variables change with distance from the coast. Taken together, these climate and geographic indicators can be interpreted to reflect the strong increase in continentality with increasing distance from the coast, which achieve maximum levels in the lee of the highest peaks of the Southern Alps. In the most interior part of the South Island, winter temperatures are lowest, summer has the highest potential evapotranspiration rates and temperature seasonality is greatest.  proxies of continentality, including temperature annual range (most extreme shifts in interior), December PET (with drier summers in the center of the South Island), and July T min (coldest winters in the interior). We therefore interpret the river δ 18 O values as most strongly responding to indicators of continentality on the South Island, and this interpretation is confirmed by the strong correlation between river δ 18 O values and distance from the coast (r = −0.72).
Compared to physiographic variables, there are weak correlations between river δ 18 O and latitude (r = 0.28, p < 0.05) and longitude (r = 0.16, p < 0.05), demonstrating that geographic variables are less important than temperature and continentality controls on δ 18 O values (Table S3). We attribute this finding partly to the NE-SW orientation of the South Island and the Southern Alps, but also to the lack of a significant latitude effect as witnessed by similar values along the entirety of the Tasman Sea coast. The correlation to sample collection altitude (r = −0.55) was weaker than those related to the temperature variables and distance from the coast, but stronger than precipitation variables. However, collection altitudes are biased to the altitudes of the road network and are thus not likely a robust proxy for the mean catchment altitude, which could have a different correlation to δ 18 O values. To overcome this bias, we also determined maximum and upslope mean hypsometric catchment altitudes for a subset (n = 151) of the data, excluding the larger rivers, for which ArcGIS could accurately identify drainage basins. Some of the small streams in our data set did not result in accurate catchments due to the relatively coarse resolution of the elevation model. The correlations between δ 18 O and mean and maximum catchment altitudes are r = −0.41 and r = −0.28, respectively, which are surprisingly weaker than to those of sample collection altitude (r = −0.51 for the n = 157 subset). The correlation between δ 18 O and catchment relief (maximum catchment altitude minus sample collection altitude) was r = −0.14, a relatively weak relationship.
We evaluated the spatial distribution of deuterium excess (d-excess) values obtained for the Tasman and Pacific slope drainages separately to determine if moisture source is another possible climatic control on river isotope values. This separation by drainage aspect was suggested by the observation of higher d-excess values (±one standard deviation) of rivers draining the Tasman (d-excess = 11.0 ± 2.6‰) versus the Pacific coast (d-excess = 8.2 ± 2.9‰; Figure 5). The origin of the different d-excess signature on the two slopes is somewhat unclear but is likely related to precipitation processes, as the strongest predictors for d-excess (Table S3) among the climatic variables were precipitation of the coldest quarter (r = 0.48), the aridity index (r = 0.42), and June precipitation (r = 0.40). Higher PET in the lee of the Southern Alps may also produce more evaporation of river values and result in higher d-excess values. Temperature variables were weak predictors of river d-excess (r ≤ |0.23|). The differences in d-excess as a function of Tasman versus Pacific slopes thus seem to indicate differences in moisture sources (Pfahl & Sodemann, 2014).
We followed this analysis by determining Tasman versus Pacific slope surface water lines and altitude effects, as well as investigating whether the difference in isotopic values was related to the presence of glaciers in the river headwaters. The respective water lines for the Tasman and Pacific slope surface water line subsets are δ 2 H = 7.18 × δ 18 O + 5.7 and δ 2 H = 7.67 × δ 18 O + 5.4, respectively. The regression equation slopes >7 indicate that a kinetic isotope effect from evaporation is not a strong control on the d-excess differences. In comparison to the relatively weak correlation between δ 18 O and sample collection altitude for the entire data set (r = −0.51), we observe stronger altitude-δ 18 O correlations when the Tasman and Pacific slope data LACHNIET ET AL.
10.1029/2021PA004220 9 of 20 are treated separately: Tasman slope streams of Δδ 18 O = −2.17‰ per km (r = −0.66; Figure 6), and the altitude effect for the Pacific slope drainages of Δδ 18 O = −2.23‰ per km is similar but has a weaker correlation (r = −0.50). These data support the contention of an altitude influence on δ 18 O values, but this relationship is obscured somewhat in the full data set because of the decrease in δ 18 O across the South Island associated with rainout and moisture source mixing. We do not find any clear difference in glacierized versus nonglacierized catchments with respect to δ 18 O and d-excess (not shown), suggesting that the presence or absence of glaciers is not a strong control on the difference in isotopic values of the Tasman and Pacific slopes.
Finally, we compare the spatial distribution of d-excess values to the extent and boundaries of South Island regional climate districts (Figure 7) to evaluate potential mechanisms that drive the observed distribution in surface water isotopes (Kidson, 2000;Lorrey et al., 2007). Salinger and Mullan (1999) applied rotated principal component (RPC) analysis to New Zealand temperature and precipitation station data from the 1930s to the 1990s to identify characteristic patterns of interannual variability. RPCs were then correlated to regional atmospheric circulation indices to evaluate how the interaction of atmospheric flow with topography produces characteristic climate zones. Four RPCs with significant correlations to atmospheric circulation indices were identified for the South Island. Two of them show how positive precipitation anomalies along the eastern margin of the South Island occur due to strong easterly to southeasterly flow derived from the Pacific Ocean. A third RPC incorporates most of the South Island west of the Southern Alps, where anomalously high precipitation is derived from westerly to northwesterly flow over southern New Zealand. The final RPC showed significant correlations between rainfall over the northernmost margin of the South Island and anomalous northerly flow. Kidson (2000) expanded on that work by identifying twice-daily synoptic-scale circulation patterns defined by k-means clustering, and retained three climate zones for the South Island (Zones 4-6; Northern, Eastern, and Western/southern South Island; NSI, ESI, and WSI). These homogenous regional climate districts provide an interpretive framework for evaluating the d-excess variations in South Island surface waters because they are associated with distinct moisture source regions. Low d-excess values (7.0 ± 2.2) along the eastern margin of the South Island are likely caused, at least in part, by precipitation derived from the Pacific Ocean along an easterly flow trajectory. Higher d-excess values (10.3 ± 4.2) west of the Southern Alps in the WSI are likely due to precipitation that originated south and west of NZ. For the northwestern NSI, d-excess values (10.8 ± 2.0) are statistically identical to the WSI, an observation that suggests that precipitation at northernmost margin of the South Island has a northerly Tasman Sea source. We argue that different moisture source conditions, specifically differences in temperature and relative humidity where surface ocean water evaporation takes place (Pfahl & Sodemann, 2014), contribute to the variations in d-excess around the South Island. Small changes in d-excess could also arise from evaporation in rivers or soil between the source regions, but this effect would have to be small because the slope and intercept of the surface water line (of δ 2 H = 8.17 × δ 18 O + 10.57) indicates minimal evaporation of waters.

Climatic and Topographic Control of South Island River δ 18 O Values
To create a predictive model of South Island river δ 18 O values as measured during the November 2016 sampling campaign, we tested several multiple regression models for δ 18 O and evaluated their statistical fits using the root-mean-square error (RMSE) metric. We screened the data set for only those correlates LACHNIET ET AL.
10.1029/2021PA004220 that explained more than 5% of the variation. We limited our model to the first four strongest variables to avoid over-fitting and collinearity, but inclusion of additional variables would yield marginally stronger correlation statistics. To provide a robust model RMSE and regression coefficients, we completed 1000 multiple regression iterations in a jackknifing approach by withholding randomly selected subsets comprising 20% (n = 38) of the data set.
First, we completed a multiple regression on the subset of samples (n = 151) that included the catchment parameters of hypsometric mean, maximum, and minimum altitudes, and catchment relief (maximum minus minimum altitude above sample location). This jacknifed multiple regression included the variables April vapor pressure, annual evapo-transpiration, mean annual temperature, and catchment relief and returned r = 0.89 (p << 0.05). Because the creation of an input grid for isoscape construction is computationally expensive (hypsometries for ∼1.5 million pixels would need to be generated to create an input grid for isoscape construction), the catchment parameters were not considered further for the isoscape model at this time; future work could explore this idea in more detail. Second, we used variables with suitable input grids in a jacknifed multiple regression that included four variables (November solar radiation, December potential evapo-transpiration, minimum temperature of the coldest month, and latitude), which has similarly high explanatory power (r = 0.88, p << 0.05) and accounts for 77% of the variance in South Island river δ 18 O values: The RMSE based on the withheld data is 0.68‰. Bivariate plots of the input variable correlations and the final model correlation to measured δ 18 O values are shown in Figure 8. The multiple regression model captures the ∼8‰ variation across South Island river δ 18 O values. A map of the residuals shows that sites with residuals between +1 to 2‰ are evenly distributed across the South Island, whereas there is a general clustering of sites with residuals between −1.0 and −2.0‰ around the highest portion of the Southern Alps ( Figure S3). We suspect that the larger number of samples with residuals that exceed ±1‰ around the Southern Alps may reflect the more complicated hydrology associated with snow, glacier ice, and liquid precipitation; however, there is no clear relationship between mean residual values or distributions for glacierized versus nonglacierized basins, suggesting that the cause is not due to the presence or absence of glacial ice in the basin catchments. In our model, the residuals were uncorrelated to the modeled δ 18 O variations ( Figure S4a), suggesting that the model is robust to the predicted δ 18 O variations. However, because our multiple regression model contains residuals that are correlated to the measured δ 18 O values ( Figure S4b), we applied a "hybrid" approach (Bowen, 2010) to produce a residual-corrected isoscape of the November 2016 river δ 18 O values (Figure 9). To create the residual-corrected map, we linearly interpolated the residuals across the South Island and tapered them to zero along the coastline, then added back the residuals. Both the original and the residual-corrected isoscapes display the same geographic variations, so our following interpretations are not affected by the choice of residuals treatment; both models are available in GeoTIFF format as supplementary information ( Figure S5), as is a derivative δ 2 H isoscape ( Figure S6).

Model Validation
The utility of any model is assessed by its ability to predict river δ 18 O values collected in different places or different times. Because there are many possible influences on river δ 18 O values, it is difficult to test the model utility based solely on the internal statistics. To illustrate, we consider that if there are significant local watershed effects including variable runoff, ages, or water sources (snow vs. rain vs. groundwater), or large seasonal variability in the river δ 18 O values, and so on, our model may not have wide applicability beyond the November 2016 sample interval and we would expect to have large differences in predicted versus measure stable isotopic values. This case would represent the "signal" being smaller than the local scale "noise." On the other hand, if our model can predict river isotope values within a reasonable uncertainty then it would suggest that our data have wider applicability beyond the interval of the sample campaign, LACHNIET ET AL.   and that the dominant climatic and physiographic controls on our data set are applicable on longer time scales, that is, the signal is stronger than the noise. With this question in mind, we ask how well is our model-based on samples collected over a short time interval of November of 2016-able to predict the δ 18 O values of independently collected archived samples across all seasons spanning the 1980s through 2010s? We first look to a subset of rivers (n = 24) collected from the South Island (Table S2) that were sampled 10 times over a year, to assess the degree of seasonal δ 18 O variation. The standard deviations of the seasonally varying δ 18 O values for the individual rivers range from 0.1 to 1.6‰, with a mean standard deviation of 0.33‰. Taken together with the 0.3-0.6‰ seasonal standard deviations of the Stewart et al. (1983) river data, these data show that the expected seasonal variations of δ 18 O values in New Zealand streams should be less than ∼1.0‰ for most rivers, and considerably less for some others. Because the signal in our spatial δ 18 O variations of ∼8‰ is sufficiently larger than ∼1.0‰ possible seasonality, then the relative spatial variations are likely robust.
To further assess model performance, we compiled 127 independently measured streams with δ 18 O and/or δ 2 H values from the published literature (Blackstock, 2011;Kerr et al., 2015;Stewart et al., 1983) including the 24 rivers reported above (Table S2) to calculate the statistical fit, and mean and standard deviation of the residuals. The Stewart et al. (1983) δ 2 H data, collected in the 1970s and 1980s, were converted to δ 18 O values according to our local surface water line equation from Figure 2. We omitted the Waimakariri River data from Blackstock (2011) because it is a large river sourced in the Alps that was sampled near the Pacific coast, but we used the other smaller rivers in this data set. The river samples were collected at different months within the year, and thus contain variability associated with the seasonal cycle, and their measured versus predicted values provide a preliminary estimate of how well our statistical model can predict the δ 18 O values of South Island rivers. We emphasize that our model and isoscape is publicly available to permit additional testing against archived or future river isotope values by other workers, which would be particularly helpful should currently unpublished data eventually be made publicly accessible. First, we compared the measured δ 18 O values of 30 rivers analyzed by Stewart et al. (1983) and Kerr et al. (2015) for which we also LACHNIET ET AL.

10.1029/2021PA004220
13 of 20 have measured δ 18 O data. The mean δ 18 O difference for this comparison is 0.09 ± 0.68‰. This comparison suggests that there have been only small variations in river isotope values over three to four decades, at least for this data subset. The small δ 18 O differences support our observations (above) that seasonal δ 18 O variation in South Island streams were also small.
Second, we interpolated model δ 18 O values for each sample in the validation data set using the δ 18 O isoscape ( Figure S4) derived from the residual-corrected multiple regression model, and calculated residual δ 18 O values for every river in the validation data set in Figure 10   Finally, we compare the river δ 18 O isoscape to meteoric precipitation (Figure 9). This exercise is not to suggest that river water sampling campaigns should be used to infer meteoric precipitation δ 18 O-that should be done by precipitation sampling-but it serves to test how well the robust spatial variations recorded by rivers are related to patterns in meteoric precipitation. As discussed earlier, if the rivers are strongly affected by seasonality and basin-scale complications, their δ 18 O values should bear little resemblance to meteoric precipitation δ 18 O values. However, if the river values encode primary climatic information that is robust over time, they should compare well to meteoric precipitation δ 18 O values. Thus, we compared the δ 18 O values of annual and seasonal means meteoric precipitation from the CDRP data set of Frew et al. (2011) to our river δ 18 O multiple regression model, by interpolating predicted δ 18 O surface water values for the locations of the isotope monitoring stations that have data spanning 2007-2010 (Table S4). The RMSE of the predicted values from the surface water isoscape versus the arithmetic mean of precipitation δ 18 O for the South Island stations is 0.68‰ (r = 0.83, p << 0.05; Figure 10), but that for weighted precipitation was considerably worse (RMSE = 1.49, r = 0.71, p << 0.05). We then compared the arithmetic and weighted mean meteoric precipitation δ 18 O values of June through November (JJASON) to the surface water isoscape values interpolated to the meteoric station locations, under the supposition that these months approximate the months that are most relevant for the November sampling dates of rivers. The statistics showed a good fit, with RMSE values of 0.79 and 0.76 for the arithmetic and weighted JJASON means, respectively (r = 0.82 and 0.85, respectively, both p << 0.05).
This model validation exercise demonstrates that our South Island river water multiple regression model and isoscape can predict δ 18 O values of independently collected rivers to better than 0.80‰, and the river data are strongly related to the seasonal (JJASON) meteoric precipitation δ 18 O values to better than ± 0.80‰. For comparison, the island-wide δ 18 O variability we see in our data set span ∼8‰. These data conclusively show that the spatial river δ 18 O values in our data set contain strong climatic and physiographic controls that are also present at other times and places in surface and meteoric waters.

Discussion
Our study demonstrates that South Island river δ 18 O values are controlled by climate and topographic effects, and these data can help constrain the dominant climatic processes that influence δ 18 O in modern and ancient settings. The ∼7‰ magnitude of δ 18 O decrease in the isotopic rain shadow is evident by comparison of δ 18 O values from small drainages on the Tasman slope of ∼−5‰ (e.g., the Punakaiki, Nile, and Waimangaroa rivers) to ∼−11 to −12‰ in the isotopic rain shadow (e.g., Motatupo, Sonora Creek, and Luggate Creek). We suggest that the large temperature drop from the Tasman sea to the region in the lee of the Southern Alps drives isotopic distillation of heavy isotopes out of air masses, which in turn produces the lowest values in areas with the lowest minimum winter temperatures. Furthermore, the region with coldest winters in the central South Island (Central Otago) is also the least affected by delivery of moisture from the Pacific Ocean with higher δ 18 O values.
There are strong winter temperature correlations with river δ 18 O values, with lowest δ 18 O values in the most continental region with the coldest temperatures. The δ 18 O-temperature gradients of 0.58-0.65‰ /°C (for July and August minimum temperature, respectively) are similar to the gradient for global stations in the 0 °C-20 °C temperature range of 0.58‰ /°C (Rozanski et al., 1993). Similarly, the isotopic temperature gradients closely bracket the gradient from midlatitude sites in winter meteoric precipitation, of 0.58‰ /°C (Fricke & O'Neill, 1999). The similarity in the isotopic temperature gradients for the midlatitude New Zealand rivers and for the midlatitude meteoric precipitation is strong evidence of a winter temperature control on South Island river δ 18 O values.
The river data expand upon and complement early work on precipitation and surface water samples that demonstrated the existence of an isotope rain shadow in the lee of the Southern Alps (Chamberlain et al., 1999;Kerr et al., 2015). We measured river δ 2 H in the isotopic rain shadow as low as −80 to −90‰, significantly lower than the −72% measured by Chamberlain et al. (1999). We have shown that the development of the isotopic rain shadow is associated with indicators of continentality, which include high temperature seasonality, low rainfall, and increased distance from the coasts. These continentality indicators have their strongest expression for all of New Zealand in central Otago and Canterbury in the lee of the Southern Alps (Macara, 2015). This situation arises from the rain shadow effect associated with rainout of westerly air masses and the foehn effect that develops from descending air that warms as it passes over the Southern Alps main divide and travels downslope. The driest valleys in the isotopic rain shadow have a semi-arid climate with median annual precipitation amounts less than 600 mm/yr, compared to >6,000 mm/yr along the crest of the Southern Alps (Macara, 2015). Our map of δ 18 O variations (Figure 1) shares similarities with the precipitation isoscape derived from a multiple regression of precipitation-weighted mean annual isotope values against temperature, precipitation amount, and elevation (Baisden et al., 2016). Both maps capture heavy isotope depletion in the lee of the central Southern Alps and higher values along the Tasman and Pacific coasts. However, differences in sample density across the South Island results in more detailed spatial variability in the surface water data.
The wide spatial extent of the isotopic rain shadow suggests that shielding of inland central Otago by coastal mountain ranges limits the delivery of easterly Pacific-sourced moisture, and that the source of river waters in the rain shadow is dominated by the Tasman Sea. A smaller area of heavy isotope depletion is predicted in the northeast South Island in the Kaikoura region where several topographic peaks exceed 2,000 m in altitude ( Figure 9). Higher δ 18 O values along the eastern coasts of the South Island indicate that an easterly Pacific moisture source is contributing less-fractionated moisture due to shorter transport distance than inland locations, reiterating previous observations (Stewart et al., 1983).
Because of the strong correlations between winter temperature and river δ 18 O, paleoclimate proxy records that are linked to South Island surface waters, such as open-system lakes with hydrogen-bearing lake sediment (for δ 2 H analyses) or carbonates (for δ 18 O) should thus be sensitive to winter minimum temperatures. Aquatic and riparian vegetation that taps into river water should also be sensitive to temperature-driven variations in δ 2 H and δ 18 O values. Similarly, sites with the largest degree of temperature seasonality and temperature annual range are also inversely correlated with lower river δ 18 O values. The strong correlations to winter temperatures suggest that isotopic proxy records associated with river systems may be robust paleotemperature indicators.
Our data have implications for the site suitability of paleoclimate proxy data. Paleoclimatic records from the South Island have been used to reconstruct climatic variability and atmospheric regimes related to increased frequency of low pressure systems embedded in this prevailing westerly flow (Lorrey et al., 2014(Lorrey et al., , 2007. However, the sensitivity of different locations on the South Island to variability of the westerlies is underdefined, and generation of δ 18 O and δ 2 H isoscapes may be used to better guide site selection for modern and paleoclimate isotopic studies and provide a means to make more informed paleoclimate interpretations from downcore data. The δ 18 O data indicate that sites along the east coast reflect a mix of low-δ 18 O Tasman-and high-δ 18 O Pacific-sourced moisture. Changes in the relative amounts of westerly versus easterly sourced precipitation over time could thus have a strong influence on δ 18 O-based climate and environmental change proxies located east of the Southern Alps. Such changes in the atmospheric flow have been noted in the observational record (Salinger & Mullan, 1999), where easterly wind anomalies are associated with higher rainfall east of the Southern Alps. Rainfall from this direction would tend to diminish the magnitude of the isotopic rain shadow because Pacific-sourced moisture experiences less orographic rainout along a trajectory into the interior of the Island. The isoscape also identifies the location having the least amount of Pacific-sourced moisture in central Otago (lowest δ 18 O values), and proxy records from the eastern margin of the δ 18 O-low region may be most sensitive to small changes in the contribution of Pacific moisture over time. Further, the amount of heavy isotope depletion may have varied over time with changes in the specific humidity of Tasman-sourced air masses. For example, speleothem isotopic evidence was interpreted to indicate that the Last Glacial Maximum (LGM, ∼21,000 years ago) had weaker winter westerlies and higher δ 18 O values than at present (Lorrey et al., 2020;Whittaker et al., 2011). Less Rayleigh distillation of heavy isotopes associated with weaker westerlies may have resulted in a decrease in the trans-alpine isotopic gradient during the LGM.
Targeting contemporaneous paleoclimate proxy sites from the Tasman slope and isotopic rain shadow regions east of the main divide may help further constrain changes in the spatial δ 18 O gradient over time and allow better testing about hypotheses of westerly wind strength and changing moisture sources. Comparisons to reactive transport models that explain δ 18 O gradients as functions of topography and energy balance appear well suited to quantify the possible causes of spatial changes in δ 18 O at different times (Kukla et al., 2019). For example, changes in effective precipitation recorded at a network of proxy sites located on different sides of the main divide, and especially eastern regions which incorporate signals of wetting and drying relative to present-day conditions reflecting changes in continentality, can support detailed assertions about synoptic type and climate regime changes (e.g., Lorrey et al., 2007Lorrey et al., , 2008. For example, an increase in blocking regime occurrence (Kidson, 2000), which is known to produce more frequent northerly and easterly flow across the country, may change the spatial isotopic gradient at times of different circulation dynamics. The minimum δ 18 O values in the isotopic rain shadow might also change over time as moisture sources shift, for example, a more easterly moisture source delivering precipitation to inland areas would plausibly increase the δ 18 O values of surface waters because there are less orographic barriers in the transport pathway compared to air masses that pass over the Southern Alps. For records that can preserve deuterium excess variation like ice cores, such data could also implicate changing proportions of Tasman Sea versus Pacific Ocean sources over time.
Our focus in this paper was primarily on the relative spatial variations in river δ 18 O values, and in the quantification of δ 18 O-temperature gradients. Our data represent a spatial "snapshot" of river δ 18 O values that can be used as a baseline to test for changes over time with future sampling campaigns. Evidence from previous measurements on river δ 18 O and comparison to our data suggests that absolute δ 18 O values are accurate to generally better than ±1.0‰, but the spatial patterns-in particular, the low δ 18 O values in the Southern Alps rain shadow-are likely to be robust to small changes in seasonality and/or climate change. Because we currently lack comprehensive island-wide data constraining the seasonal or inter-annual differences in river water δ 18 O values, the isoscape patterns should be tested and refined against future sampling campaigns which may improve our understanding of isotope climatology processes. Such a comparison would be useful, for example, to test for changing moisture sources under regional warming associated with anthropogenic climate change. In addition, strong spatial signatures may arise due to changes in major ocean-atmosphere modes of variability that impact on New Zealand rainfall like the El Niño/Southern Oscillation (Jiang et al., 2013), the Southern Annular Mode (Ummenhofer & England, 2007) or just interannual δ 18 O variability. The δ 18 O data presented in this work are best viewed as a first attempt to determine river water spatial variations upon which future work can build. Our comparison with limited isotopic data from the 1970s and early 1980s suggests relatively small changes in isotopic values over time. However, a more robust test of anthropogenic climate influenced on river water isotope values would require resampling of the network over time. We hypothesize that δ 18 O values would increase over time due to the temperature-dependent fractionation effect as the South Island warms due to a poleward movement of westerly storm tracks and expansion of the subtropical ridge, as is predicted for future drying of New Zealand due to a more positive Southern Annular Mode under the high-emissions RCP8.5 scenario (Lim et al., 2016;Yin, 2005).
Finally, our data also show the utility of short duration but areally extensive surface water sampling campaigns to constrain climatic controls on oxygen and hydrogen stable isotopic values. Despite the limited time interval for such sampling campaigns, such data reveal robust spatial variations and that have predictive ability. The conclusion that spatial surveys of surface water isotopes encode primary climatic controls is remarkable given the wide variety of potentially confounding processes that may influence surface water δ 18 O values. That a simple multiple regression model, collected over a short time interval, can predict δ 18 O values of independently collected samples at different times and places with a precision of better than ±0.80‰ (out of an island-wide range of 8.0‰) is strong substantiation for the validity of our approach. The data presented here for the South Island of New Zealand extend and confirm surface water tours in other regions spanning the tropics to the high-latitudes that have also shown strong predictive ability (Kendall & Coplen, 2001;Lachniet & Patterson, 2006, 2009Lachniet et al., 2007Lachniet et al., , 2016Li & Garzione, 2017). Furthermore, the surface water isoscape shows a reasonable ability to predict meteoric precipitation δ 18 O values (with an RMSE = 0.76 for weighted JJASON), an observation that is robust despite the complications that arise from river δ 18 O seasonality, residence times, and watershed hydrologic processes.