A method for predicting hydrogen and oxygen isotope distributions across a region's river network using reach-scale environmental attributes

Stable isotope ratios (isotope values) of surface water reflect hydrological pathways, mixing processes, and atmospheric exchange within catchments. Development of maps of surface water isotope values (isoscapes) is limited by methods to interpolate point measures across river networks. Catchment attributes that alter surface water isotope values affect downstream river reaches via flow, but some attributes such as artificial dams are no more likely to affect nearby unconnected catchments than distant ones. Hence, simple distance-based geospatial and statistical interpolation methods used to develop isoscapes for precipitation and terrestrial systems are less appropriate for river networks. We used a water-balance-based method to map long-term average δ²H and δ¹⁸O for New Zealand rivers, incorporating corrections using catchment environmental predictors. Inputs to the model are national rainfall precipitation isoscapes, a digital elevation layer, a national river water isotope monitoring dataset (3 years of monthly sampling at 58 sites), and river environmental databases covering around 600 000 reaches and over 400 000 km of rivers. Much of the spatial variability in δ²H and δ¹⁸O of New Zealand river water was explained using the initial combination of precipitation isoscapes and a simple water balance model. δ²H and δ¹⁸O isoscapes produced by subsequently applying residuals from the water balance model as a correction factor across the river network using regression kriging showed improved fits to the validation data compared to the correction using ordinary kriging. Predictors of high importance in the regression included upstream lake and wetland area, which was not strongly spatially autocorrelated nationally. Hence, additional hydrological process information such as evaporation effects can be incorporated into river isoscapes using regression kriging of residuals. The resulting isoscapes have potential applications in ecological, hydrological, and provenance studies that consider differences between surface water isotope values and those of other components of the hydrological cycle (e.g. subsurface runoff or local precipitation).

1 Introduction

Stable hydrogen and oxygen isotope measurements in surface water provide powerful tools for hydrology, ecology, and food science; uses include identifying runoff sources to rivers, determining migratory pathways of animals, and tracing the provenance of food (Kelly et al., 2005; Cable et al., 2011; Gao and Beamish, 1999). Typically, the utility of tracers in mixing and source partitioning studies relies on the tracer concentration in potential sources being sufficiently distinct, relative to the degree of statistical noise in the system (Fry, 2006). Therefore, it is important to know the spatially explicit distributions of hydrogen and oxygen stable isotope ratios (hereafter isotope values following delta notation; Coplen, 2011) for waterbodies. Isotope values of river water are strongly dependent on those of precipitation (Kendall and Coplen, 2001), mixing of water sources, and isotopic fractionation resulting from different hydrologic processes occurring at and below the land surface within their catchments. Variations in the isotope values of precipitation reflect the temperature at which condensation occurs and prior condensation as air masses are transported over land, resulting in the following:

1. Precipitation becomes more depleted in heavy isotopes at higher latitudes and higher elevations (Dansgaard, 1964).

2. Precipitation inland tends to be more isotopically depleted than that falling at the coast (Dansgaard, 1964; Winnick et al., 2014).

3. At mid- to high latitudes, precipitation isotope values vary seasonally with changing temperatures, becoming more enriched in summer and depleted in winter (Craig, 1961).

The isotope composition of precipitation at a site is also affected by the source and trajectory of water vapour in air masses (Jouzel et al., 2013; Crawford et al., 2013; McDonnell, 1988) and sub-cloud processes including re-evaporation (Risi et al., 2008; Winnick et al., 2014). Understanding of the processes controlling precipitation isotope values has aided development of precipitation and surface water isotope maps (Bowen and Revenaugh, 2003; Bowen et al., 2011), with resulting hydrological and cross-disciplinary applications (Vander Zanden et al., 2016; Jasechko et al., 2013).

A range of hydrologic processes occur at and below the land surface that alter the isotope composition of surface waters relative to the isotope composition of precipitation. These include fractionation during atmospheric exchange processes at ground, water, and vegetation surfaces (Evaristo et al., 2015; Ehleringer and Dawson, 1992; Gat, 1996), temporal and spatial patterns of sub-surface flow (Gonfiantini et al., 1998), and mixing between isotopically distinct waters (Klaus and McDonnell, 2013; McDonnell et al., 1991). A key limitation of the current approaches to modelling stable isotope values of river water is the representation of these processes across river networks (Bowen et al., 2011; Brennan et al., 2016). Geostatistical interpolation methods and linear regression used for the spatial modelling of precipitation are less appropriate for surface waters because some landscape attributes that alter surface water isotope values, such as dams and lakes that increase atmospheric exchange (Gibson et al., 2016; Craig, 1961), are transferred to reaches downstream but have no impact on nearby, unconnected river networks. Hence, isotope values of surface waters are controlled by a mix of geospatial dependencies and those dendritic dependencies unique to river networks (Brennan et al., 2016).

The water balance model approach of Bowen et al. (2011) provides a method of mapping hydrogen and oxygen isotope ratios across river networks that incorporates variation that can be modelled with interpolation approaches and variation that follows the dendritic patterns unique to river networks. This method takes advantage of spatial and seasonal predictability of precipitation isotope values, evaporative water loss, and catchment discharge. Spatial patterns of residuals between predicted (modelled) and measured (annual averages of measured stream water isotope values) results are then interpolated using ordinary kriging to improve models of river water isotopes. Bowen et al. (2011) note that residual adjustment using ordinary kriging only indirectly considers many of the landscape attributes that may affect stream water isotope values. These attributes include those that affect the temporal and spatial patterns of sub-surface flow, such as catchment slope (Salama et al., 1993), and those that affect atmospheric exchange processes, such as the presence of dams and natural lakes (Halder et al., 2015; Gat, 1996). When spatial variability within the catchment is low and sampling density is high, ordinary kriging methods may be appropriate; however, New Zealand is an example of a country with relatively large distances between long-term river monitoring sites but high variation in relief and climate over relatively short distances. Under these conditions, using regression-based statistical methods with databases of explanatory environmental characteristics, rather than simple distance-based methods, to explain hydrological variation across the reaches of river networks provides a superior prediction of spatial variability (Snelder and Biggs, 2002; Hicks et al., 2011). Regression kriging is a spatial prediction technique often used in soil mapping that combines spatial similarity with non-location predictors (Hengl et al., 2007; Keskin and Grunwald, 2018). This technique would seem appropriate for correcting surface water isoscapes, for which the divergence of isotope values away from those predicted using a water balance approach can be explained by catchment characteristics that are not strongly spatially autocorrelated. In this paper, we follow the spatiotemporal water balance method of Bowen et al. (2011) to derive long-term average water isotope values from New Zealand rivers, using precipitation isoscape coefficients of Baisden et al. (2016), gridded meteorological data from National Institute of Water and Atmospheric Research (NIWA)'s virtual climate station network (Tait et al., 2006; Tait, 2008), and a digital elevation map (DEM)-based river network (Snelder and Biggs, 2002). We calculate residuals between modelled surface water isotope values and values measured at 58 long-term river water sampling sites across New Zealand. We then extrapolate these residuals across New Zealand using a regression method that incorporates explanatory environmental variables from the River Environment Classification (REC) database (Snelder and Biggs, 2002), Freshwater Ecosystems of New Zealand (FENZ) geodatabase (Leathwick et al., 2010), and geographic information system (GIS) lake data as predictors. We compare this regression kriging residual estimation with ordinary kriging (Cressie, 1992; Bowen et al., 2011). We employ the method above to develop the first complete isoscape for the rivers and streams of New Zealand (NZ).

2 Methods

2.1 Water balance model for estimation of river water stable isotope values

Our initial method to calculate long-term isotope values of river water closely followed that given by Bowen et al. (2011) to calculate their “spatiotemporally weighted water balance model” and is only briefly summarised here. For each 5 km × 5 km grid cell, we estimated monthly runoff (Q) as the larger value of precipitation minus evaporation (P−E) or 0.01×P. The latter estimate (0.01×P) provides an estimate of dry season and dry region runoff (Bowen et al., 2011). We estimated the annual isotopic flux of runoff from monthly weighted values for each cell, following Eq. (1):

$$\begin{array}{}\text{(1)}& \mathit{\delta }{Q}_{\mathrm{m}}=\sum _{i=\mathrm{1}}^{\mathrm{12}}\left[\left(P_i-E_i\right)\times {\mathit{\delta }}_i\right],\end{array}$$

where δ_i is the isotope value of precipitation within the grid cell at month i. The monthly weighted isotopic flux (δQ_m) is calculated by summing across all monthly values (i). Routing of the gridded runoff started with a 30 m gridded DEM of NZ. We pit-filled the grid using Spatial Analyst in ArcGIS 10.6 and then calculated flow direction next for each grid cell using the d8 method (O'Callaghan and Mark, 1984), which considers flow to one of eight adjacent cells. Values of Q and δQ_m in each grid cell were then accumulated downstream using the Flow Accumulation tool in ArcGIS.

2.2 Model input data

Gridded precipitation isotope values for New Zealand were generated using the monthly models of Baisden et al. (2016). We applied the regression coefficients of Baisden et al. (2016) to the geographical and climate predictor data given in their study, including elevation, latitude, pressure, solar radiation, air temperature, vapour pressure, wind speed, soil temperature at 100 mm depth, and an estimate of potential evapotranspiration (PET). Daily climate variables, PET, and precipitation were retrieved from the NIWA Virtual Climate Station Network (VCSN; Tait et al., 2006; Tait and Woods, 2007), covering the period 1997–2018. Precipitation-weighted means were produced for climate drivers, and coefficients were applied to yield monthly δ²H and δ¹⁸O values for precipitation. While Baisden et al. (2016) also used VCSN interpolated data in their regressions, the VCSN is regularly updated as new physical climate stations are added to or removed from the national climate network and as the spline models are periodically improved (Tait et al., 2006). Hence, we checked the results of our procedure by performing regressions between our modelled, amount-weighted monthly precipitation isotope values and measured values from the dataset of Baisden et al. (2016), which comprises monthly collections from 51 sites across New Zealand between 2007 and 2010. Actual evapotranspiration (AET), as an input to the water balance model, was estimated according to the model of Porteous et al. (1994), assuming a soil water storage capacity of 150 mm and a critical value of 50 % of this capacity, above which AET was assumed to be equal to PET and below which PET was scaled linearly to zero. Digital elevation input data were sampled from a 30 m digital elevation model of New Zealand at the respective VCSN locations. Catchment environmental characteristics were taken from the REC (Snelder and Biggs, 2002) and FENZ geodatabases (Leathwick et al., 2010).

A summary of contributing data sources is shown in Fig. 1. Table 1 gives a summary of water isotope datasets used in model calibration and validation.

Figure 1 GIS layers and data sources used in model development. Panels (a)–(f) show GIS layers used for water balance modelling, plotted over the Canterbury region in New Zealand's South Island. Precipitation, actual evapotranspiration (AET), and stable isotope values of precipitation (a–c) are used to calculate catchment-scale, flow-weighted isotope values of recharge (δQ; i.e. flow × isotope δ value). The river network (f) is derived from the DEM (d) and flow direction (e) layers by setting a suitable flow accumulation threshold. Catchment δQ values are routed downstream through the river network to generate reach-scale estimates of δ²H and δ¹⁸O. Panel (g) shows river water isotope measurement sites (triangles) used for model correction and major lakes (blue area).

Table 1 Water isotope datasets used in model calibration and validation. Lake samples provided in Stewart et al. (1983) were not included in our analyses. For Stewart et al. (1983), Marttila et al. (2017), and unpublished data, annual averages for sites were used where available. Data used for river model validation are compiled in the Supplement (Sect. S1).

2.3 Calculation and kriging of model residuals

Modelled river water isotope values were compared to annual average values from 58 sites from the National River Water Quality Network (NRWQN; selected to represent catchments nationally; Yang et al., 2020). Design of the NRWQN is described by Smith and McBride (1990), while descriptions of physical (catchment), flow, and chemical conditions at monitoring sites can be found in Davies-Colley et al. (2011), Julian et al. (2017), and Yang et al. (2020). Measurements from this network have been used to develop and calibrate a range of hydrological and water quality models (e.g. Alexander et al., 2002; Elliott et al., 2005). The residual calculations used in this study compared predicted flow-weighted annual average river water isotope values to the averages of measurements. Measurements were flow-weighted monthly values at each site from April 2017–March 2020. Data from April 2017–March 2019 are published in Yang et al. (2020). For the current analysis, an additional year of data (April 2019–March 2020) was collected from all NRWQN sites following the methods of Yang et al. (2020). Briefly, samples were taken monthly from each of the 58 sites, and river flow rates were recorded at the time when samples were collected. All samples were stored in 100 mL tubes in insulated ice bins, keeping samples at approximately 0 ^∘C and in darkness in the field, and then frozen and stored at approximately −20 ^∘C in the laboratory and later thawed for analyses. Isotope analyses were conducted using an isotope ratio infrared spectroscopy (IRIS) on a wavelength-scanned, cavity ring-down spectrometer (WS-CRDS) model L1102-i (Picarro, Inc., Sunnyvale, California, USA). Any variation in isotope values resulting from diurnal patterns was minimised by the collection of samples at similar times of the day. River flow rates were recorded at the time when samples were collected.

The locations of monitoring stations were adjusted in ArcGIS to coincide with the nearest grid cell of the stream network. Modelled isotopic compositions were extracted from each of the grid cells corresponding to a monitoring site, and the differences between the modelled and observed were calculated.

To improve reach-scale predictions of water stable isotope values beyond the predictions of the water balance model, we used and compared two kriging methods, namely ordinary kriging and regression kriging. These two methods were used to extend calculated residuals between simulated and observed stable isotopes at our 58 validation sites across the New Zealand stream network. Ordinary kriging, following the method of Bowen et al. (2011), only relies on spatial similarity (i.e. spatial autocorrelation) of the residuals. Regression kriging relies not only on spatial similarity but also on non-location predictors. In this case, the non-location predictors are catchment environmental variables available at sub-catchment scale across the New Zealand from the REC and FENZ geodatabases (Snelder and Biggs, 2002; Leathwick et al., 2010).

In ordinary kriging, the prediction is a weighted mean of the following observations:

$$\begin{array}{}\text{(2)}& y_i=\sum _{k=\mathrm{1}}^n{\mathit{\beta }}_k\left(i\right)y_k^{\mathrm{o}},\end{array}$$

where y_i is the prediction at location i, $y_k^{\mathrm{o}}$ is the kth observations (k=1, …, n), and β_k(i) is the kriging weight of $y_k^{\mathrm{o}}$ to y_i, which is a function of locations of i and observation k and depends on the spatial autocorrelation structure of the variable (Hengl et al., 2007).

In regression kriging, the prediction is the sum of a simple regression (normally linear regression as shown below), and the result of ordinary kriging as follows:

$$\begin{array}{}\text{(3)}& y_i=\sum _{j=\mathrm{1}}^p{\mathit{\alpha }}_j{x}_j+\sum _{k=\mathrm{1}}^n{\mathit{\beta }}_k\left(i\right)e_k^{\mathrm{o}},\end{array}$$

where x_j is the jth environmental variable (j=1, …, p), α_j is the corresponding regression coefficient, and $e_k^{\mathrm{o}}$ is the residual at kth observation after the regression.

The application of regression kriging includes two steps. The first step is to apply a simple multivariable regression aiming to select the best independent non-location variables. An initial list of over 30 potential independent variables for the regression model were selected from the REC and FENZ geodatabases, based on our understanding of the hydrological processes which potentially influence stream isotope components, including geological factors (e.g. slope, elevation, aspect, drainage density, and lake coverage), climate factors (e.g. precipitation and evaporation), and land cover (Yang et al., 2019, 2020). Table A1 includes descriptions of these variables. Because many of these variables represent characteristics of the upstream catchment, they provide a representation of catchment-specific upstream processes (such as evaporation from wetlands, reservoirs, and lakes) that may affect river water isotope values. However, given the relatively small number of validation sites available (58), correlation analysis and stepwise regression were applied to reduce number of independent variables. From the list of independent variables in Table A1, five were selected in the regression analysis based on BIC (Bayesian information criterion), following the “one in 10 rule” (e.g. Harrell, 2015), where one predictive variable can be included for every 10 sites in the dataset. In the second step, spatial autocorrelation is considered together with the five selected variables, following Eq. (3), to give the prediction. Similarly, since nonlinear regression requires an increased number of coefficients to be estimated, we also used linear regression in Eq. (3) to avoid overfitting.

The performance of ordinary kriging and regression kriging in predicting spatial patterns of residual values was assessed across the 58 long-term river monitoring sites using k-fold cross validation (Stone, 1974). In this study, we applied LOOCV (leave one out cross-validation) and 10-fold cross-validation to assess the performance of the two kriging methods.

2.4 Data for final model validation

Final model validation used water isotope data collated from previously published studies of rivers across New Zealand (Lachniet et al., 2021; Kerr et al., 2015; Stewart et al., 1983; Marttila et al., 2017). The report of Stewart et al. (1983) includes δ²H (but not δ¹⁸O) samples taken from nearly 200 sites throughout New Zealand between 1966 and 1981. These are largely single samples for rivers but include repeated sampling over several years at some sites that highlighted storm-to-storm and seasonal variation. The work of Lachniet et al. (2021) is based on a single, high spatial resolution sampling campaign across the South Island of New Zealand in 2016. Kerr et al. (2015) reported single samples taken from small rivers across an east–west transect through the Southern Alps in the South Island. Marttila et al. (2017) conducted over 2 years of monthly sampling from seven river sites in a small area of the South Island. The literature validation dataset is provided in the Supplement (Sect. S1) and detailed in Table 1. Because this validation dataset contains largely synoptic sampling data, to increase the number of annual average values in the validation dataset, we included unpublished δ²H and δ¹⁸O data from nine river sites for which regional government staff have conducted monthly sampling over date ranges between 2017 and 2020. Water samples from these sites were stored and analysed, following the methods of Yang et al. (2020), and the annual average values were used in regression analyses. These sites are listed under Dudley et al. (unpublished data) in the Supplement. Because the dataset of Stewart et al. (1983) contained δ²H but not δ¹⁸O data, the δ²H literature validation dataset is considerably larger (418 sites) than that for δ¹⁸O (231 sites), and the δ²H literature validation dataset also contains a higher proportion of sites in the North Island because the large dataset of Lachniet et al. (2021) is restricted to South Island rivers.

3 Results and discussion

3.1 Precipitation isotope values

Our modelled estimates of long-term precipitation isotope values are a close match to those published in Baisden et al. (2016), from which we took regression coefficients. δ²H and δ¹⁸O values become increasingly enriched in the northern and low-elevation areas of New Zealand, with the sharpest spatial gradients at the western edges of New Zealand's central mountain ranges which face New Zealand's prevailing westerly–southwesterly winds and dominant western moisture source. Regressions between modelled, amount-weighted monthly δ²H values and δ²H values in monthly precipitation collections (i.e. the 2007–2010 monthly collection time series of Baisden et al., 2016) resulted in R²=0.45 and RMSE = 10.1 ‰ for monthly predictions and R²=0.79 and RMSE = 5.1 ‰ for amount-weighted average values for the full 2007–2010 time span of the field sampling. Regressions between modelled, volume-weighted monthly δ¹⁸O values, and measured monthly precipitation δ¹⁸O values resulted in R²=0.48 and RMSE = 1.2 ‰ for monthly predictions and R²=0.80 and RMSE = 0.55 ‰ for 2007–2010 average predictions. These fit statistics differ very slightly from those published by Baisden et al. (2016), showing improved monthly fits but a poorer longer-term fit for δ¹⁸O. This is likely due to updates to the climate network data used for VCSN model input and updates to the interpolation scheme, as described above.

3.2 River water stable isotopes modelled using the water balance method

Modelled, uncorrected estimates of long-term river water isotope values for each of the 58 NRWQN stations ranged from −68.66 ‰ to −32.54 ‰ (national average −49.83 ‰) for δ²H and −9.56 ‰ to −5.67 ‰ (national average −7.56 ‰) for δ¹⁸O. Average (2017–2020) measured river water isotope values at these same NRWQN sites ranged from −77.01 ‰ to −19.75 ‰ for δ²H (national average −46.38 ‰ ) and −10.63 ‰ to −2.85 ‰ (national average −7.25 ‰) for δ¹⁸O. The uncorrected model explained 72 % of the observed variability in δ²H values between sites, and 69 % of the observed variability in δ¹⁸O values between sites (Fig. 2). Notably, however, the uncorrected model showed a tendency to underpredict both δ²H and δ¹⁸O in warmer northern and low-elevation sites (i.e. at the less negative end of the range). Model underprediction for these stations is likely due, at least in part, to the evaporative enrichment of surface waters. Since these 58 NRWQN stations extend across a large latitude and elevation gradient, rivers sampled at northern and low-elevation sites flow through catchments with substantially warmer mean annual temperatures that are likely to have led to greater evaporative enrichment. Many of the sites in northern New Zealand also have large areas of lakes and wetlands in their upstream catchments (Yang et al., 2020) contributing to greater potential for evaporation.

Figure 2 Relationship between modelled and observed δ²H (a) and δ¹⁸O (b) at 58 long-term river water monitoring sites across New Zealand. Modelled isotopes were the long-term averages of simulated isotopes computed following the water balance model approach of Bowen et al. (2011). The solid line shows the regression fit, and the dashed line shows a 1:1 fit. The modelled values tend to overpredict the river water isotope at the more negative end of the range (e.g. at high-latitude cooler southern sites) and underpredict at the less negative end of the range (e.g. warmer, northern sites).

Overprediction of both δ²H and δ¹⁸O was apparent at some sites (Fig. 2). This pattern has a precedent in areas of the continental USA with high topographic relief, where the spatially weighted model of Bowen et al. (2011), which did not consider seasonal differences in P−E, appeared to underpredict the contribution of isotopically depleted runoff from high elevations to rivers. In that study, overprediction was ascribed to the use of a model that did not consider seasonal differences in evapotranspiration; this issue was largely resolved using a temporally weighted model that considered these differences. In our study, this issue was most prominent in the South Island of New Zealand, which is characterised by high topographic relief, and high-elevation recharge sources dominate the flow of the majority of large rivers. However, in our study, we use a water balance model that considers seasonal differences in P−E but does not appear to have fully resolved this pattern. Potential reasons for this are, firstly, that the precipitation model we use overestimates the isotope values of precipitation at high elevations or in leeward, rain shadow areas; of the eight sites where predicted δ¹⁸O values exceed average measured δ¹⁸O values by >1 ‰ (Fig. 2), seven are in alpine-fed rivers on the leeward east of New Zealand. Rivers and streams leeward of the Southern Alps show isotopically depleted values characteristic of spillover of orographic precipitation from the windward west of the Southern Alps (Stewart et al., 1983; Kerr et al., 2015). While the linear regression approach employed in the precipitation model we used does not directly adjust for such orographic effects, it does to some extent capture areas of low δ²H and δ¹⁸O values in eastern New Zealand (Baisden et al., 2016).

A second potential driver of model overprediction of river water isotope values is surface water/groundwater interactions. Rivers in lowland, leeward New Zealand tend to gain a significant portion of their flow from groundwater (Yang et al., 2019). Where this groundwater is dominantly derived from high-elevation recharge, gaining streams and springs tend to be depleted in stable water isotopes (Stewart et al., 2018). The (uncorrected) water balance model we employ does not consider groundwater pathways; this may result in errors particularly in areas where gains from or losses to groundwater are considerable.

Table 2 Selected environmental variables and related information and importance ranks, p values, and t statistics in the regression for δ²H and δ¹⁸O residuals. Values for δ¹⁸O residuals are in parentheses.

Figure 3 Comparison between observed residuals and spatial predictions of residuals made using two interpolation methods, i.e. ordinary kriging and regression kriging. Observed residuals are the differences between spatiotemporal model predictions of δ¹⁸O and measurements from 58 long-term river monitoring sites across New Zealand. Comparison with estimates of residuals generated using the LOOCV (a) and 10-fold cross-validation (b) suggests the improved fit of regression kriging relative to ordinary kriging.

3.3 Kriging of isotope residuals

For regression kriging, simple multivariable regression applied to δ²H and δ¹⁸O residuals, individually, gave a final list of five independent environmental variables selected on the basis of t statistics in the two linear regressions (Table 2). These variables are site elevation and average catchment elevation (SiteElev and usCatElev), average catchment slope (usAveSlope), upstream annual rainfall variability (usAnRainVar), and upstream lake and wetland area (usLWArea). The reduced regression analyses explained over 60 % of the of the variation in the isotope residuals; R² values were 0.6 and 0.66 for the δ²H residual and δ¹⁸O residual regressions, respectively. The importance of each environmental variable (second column in Table 2), ranked by t statistic, shows that the upstream lake and wetland area, upstream annual rainfall variability, and average catchment elevation are the three most sensitive variables in both δ²H and δ¹⁸O residual regressions, although the order of importance differs between the two analyses. A possible cause for the higher ranking of upstream lake and wetland area in the δ¹⁸O regression is the greater sensitivity of the ¹⁸O component of water to kinetic fractionation effects than the ²H component (Craig, 1961; Gat, 1996). Other variables in Table A1 sensitive to isotope residuals that were not included in regression included climatic variables (e.g. annual average potential evapotranspiration, annual average actual evapotranspiration, annual average precipitation, and upstream solar radiation in summer). Including these variables could improve R² values for regressions (e.g. by including annual average precipitation as a predictor, R² for δ²H, and δ¹⁸O residuals can reach 0.76 and 0.80, respectively). However, this led to overfitting and caused the spatial isotope prediction in the following steps to be unreliable. Once the above five environmental predictor variables were chosen, regression kriging was applied to estimate the spatial patterns of residuals across New Zealand. The performance of ordinary kriging, and regression kriging were then evaluated with LOOCV and 10-fold cross-validation (CV). The results of this assessment for δ¹⁸O residuals are shown in Fig. 3. Both LOOCV and 10-fold CV led to similar results, with the regression kriging simulations scattered around the 1:1 line, while those of ordinary kriging are mostly close to a horizontal line with values around 0. There is little predictive ability for ordinary kriging, as indicated by R² values below 0.1 for both LOOCV and 10-fold CV. Regression kriging achieved satisfactory results, with R² over 0.75 for both LOOCV and 10-fold CV. These results indicate that regression kriging can be reasonably used to estimate spatial distributions of isotope residuals for rivers. It is worth noting that multiple regression contributed around 0.66 and kriging contributed around 0.1 of the total R². This emphasises the importance of the multiple regression step in this method, i.e. the first part of Eq. (3).

Figure 4 Distribution of isotope residuals (i.e. difference between modelled and observed isotopes) from ordinary kriging (a, c) and regression kriging (b, d). Positive values equate to a negative adjustment of river δ²H and δ¹⁸O values from the water balance model, as seen in leeward (eastern) areas of both the North and South islands. Ordinary kriging results in regional zones of negative isotope residuals around validation sites downstream from lakes and wetlands (such as in the central North Island), while regression kriging predicts negative isotope residuals in low-elevation western areas and in all lakes and wetlands and their downstream reaches.

3.4 Spatial distributions of model residuals

Spatial interpolation of isotope residuals was applied to the entire national river network. Residuals interpolated using ordinary kriging (Fig. 4a and c) show regional patterns; the water balance model tends to overpredict δ²H and δ¹⁸O values in rain shadow areas in the east of both the North and South islands and underpredict δ²H and δ¹⁸O values in the northern and western regions of both islands. Residuals shown in Fig. 4 are calculated as the modelled minus average measured δ values for individual river reaches, so positive residuals occur where measurements were, on average, more depleted than the model estimate. Lower isotope values in precipitation and rivers leeward of large mountain ranges are driven by fractionation during orographic precipitation processes (Stern and Blisniuk, 2002; Poage and Chamberlain, 2001), and the patterns observed across the Southern Alps indicate that the dominant source of water to rivers in leeward areas is precipitation carried over the Southern Alps by wind (Kerr et al., 2015). The overprediction of δ²H and δ¹⁸O values (of the water balance model) in rain shadow areas indicates either greater recharge from high elevations to river flow within these regions or local differences between actual and predicted (modelled) precipitation isotope values. We suggest that at least some of the patterns of residuals for the water balance model (Fig. 4) can be accounted as a correction of the precipitation isotope model of Baisden et al. (2016) with regard to west–east depletion of precipitation isotope values associated with orographic rainfall across the Southern Alps. The basis of this conclusion is shown in Fig. A1, in which we have plotted residuals between measured and modelled δ²H values for small rivers across an east–west transect through the Southern Alps (Kerr et al., 2015). The precipitation isotope model (i.e. Baisden et al., 2016) appears to underpredict δ²H values to the west of the main divide and overpredict them to the east. The patterns seen in Figs. 2 and A1 indicate that these prediction errors were improved by the regression kriging correction, though not fully resolved. Spatial patterns of the strongest predictors are shown in Fig. A2; we believe that relationships between residuals and some predictors may reflect issues with representation of orographic effects by the precipitation model. All predictors in Table 2, except the upstream lake and wetland area, show a strong west–east gradient; for example, areas to the east of the Southern Alps, where the water balance model overpredicts river water δ²H and δ¹⁸O values, have lower average catchment slopes and higher upstream annual rainfall variability than that present to the west of the Southern Alps, where the water balance model underpredicts river water δ²H and δ¹⁸O values.

Figure 5 Comparison of the literature values of water δ²H across New Zealand rivers, with predicted long-term average δ²H values from the uncorrected water balance model (a), ordinary-kriging-corrected model (b), and regression-kriging-corrected model (c). Panel (d) shows a comparison between predictions of the regression-kriging-corrected model with mean annual measured values at sites with >1 year of monthly data. Lines show a 1:1 relationship.

The significance of upstream lake and wetland area in the residuals regression can, however, be attributed with some confidence to evaporation from surface waters with long durations of exposure to the atmosphere. Wetland systems and lakes are spread relatively well throughout New Zealand (Cromarty and Scott, 1996); glacier-carved lakes and artificial lakes used for hydroelectric power are common through much of the Southern Alps region of the South Island, while volcanic crater lakes are present particularly in the central–eastern regions of the North Island (Fig. 1g). Notably, the presence of validation sites on rivers draining major lakes creates areas of residuals that differ from regional trends. For example, within the central North Island region, validation sites are sparsely spread, and roughly half of the validation sites sit on major rivers draining two caldera lakes, namely Lake Taupo (616 km²) or Lake Tarawera (41 km²). Using ordinary kriging, these sites create an area of negative residuals (e.g. Fig. 4c) that are applied as a correction to all nearby river reaches. This effect suggests that the application of simple distance-weighted kriging of residuals may be inappropriate for our validation dataset, as discussed in Sect. 3.3, where sites are sparse relative to both land area, patterns of lake distribution, and topographic relief. Regression kriging appeared to account to some extent for regional-scale effects; regional patterns of overprediction in eastern rain shadow areas are represented, as are patterns of underprediction in northern and western New Zealand. However, regression kriging does not appear to create a regional pattern of negative residuals associated with lake evaporation in the central eastern North Island (Fig. 4d).

Figure 6 Predicted δ²H and δ¹⁸O of surface water of New Zealand. Panels (a) and (d) show uncorrected values from the water balance model. Panels (b) and (e) give residual-corrected values based on ordinary kriging (OK correction). Panels (c) and (f) give residual-corrected values based on regression kriging (RK correction).

3.5 Final model validation

To validate the results from the original water balance model (uncorrected), ordinary kriging (OK corrected) and regression kriging (RK corrected), we compared the model results to collated literature values from previously published studies of rivers across New Zealand (Lachniet et al., 2021; Kerr et al., 2015; Stewart et al., 1983; Marttila et al., 2017) and unpublished data from nine sites, as described in Sect. 2.4. The validation analysis required assigning a reach identification number from the REC to each record reported in the citations. Using this approach, linear regression models between RK-corrected isotope predictions and all the corresponding literature validation data explained 75.4 % and 71.3 %, respectively, of the δ²H and δ¹⁸O variance. We note that data from six sites sampled by Lachniet et al. (2021), for which we could not confidently assign a reach identification number, were not included in these regressions. The fit for δ²H is shown in Fig. 5c. Root mean squared errors (RMSEs) were 6.2 ‰ for δ²H on 416 degrees of freedom and 0.67 ‰ for δ¹⁸O on 228 degrees of freedom. Datasets incorporating repeat sampling of stable isotopes in New Zealand river water are scarce; only 23 of the 418 δ²H records in the literature validation dataset represent annual averages of monthly sampling. The river samples from the South Island rivers dataset of Lachniet et al. (2021) were collected within the month of November, within the austral spring snowmelt period; South Island river water isotope values at this time are generally below the annual average (Yang et al., 2020; Marttila et al., 2017). Hence, we would expect our model predictions at sites from high latitudes with snowmelt influence to be generally more isotopically enriched than our literature validation dataset. This effect appears strongest in the uncorrected model results. When smaller rivers (Strahler order < 3, where larger seasonal fluctuations in isotope values would be expected) were removed from the analysis, the fit of the final model improved, explaining 80.0 % and 74.2 %, respectively, of the δ²H and δ¹⁸O variance across the literature validation dataset, with RMSEs of 5.7 ‰ for δ²H on 329 degrees of freedom and 0.62 ‰ for δ¹⁸O on 175 degrees of freedom. The most appropriate sites for validation of our model were those 23 sites for which long-term monthly sampling records are available, enabling us to compare predicted and measured annual average values. Average values derived from monthly measurements made over several years are likely to average out much of the temporal variation in river water isotope values which results from temporal variation in precipitation isotope values, evaporation, and contributions from different flow pathways. At these sites, the model explained 90.6 % of the δ²H variation across the dataset. The RMSE was 2.99 ‰ for δ²H on 21 degrees of freedom (Fig. 5d).

While comparison of modelled annual average values to the full validation dataset containing predominantly single-sample values from river reaches (Fig. 5a–c) shows relatively poor model fits compared to those for long-term measurements (Fig. 5d), these comparisons serve to demonstrate an improvement from the uncorrected and ordinary-kriging-corrected models to the final regression-kriging-corrected model.

3.6 Patterns and drivers of final model performance

3.6.1 Patterns of precipitation isotopes, runoff, and evaporation from surface water are the dominant drivers of δ²H and δ¹⁸O values in rivers

Across New Zealand, spatial and temporal patterns of δ²H and δ¹⁸O in precipitation and runoff are dominant drivers of δ²H and δ¹⁸O in river water (Fig. 6; see Yang et al., 2020). The (uncorrected) water balance model, which explicitly represents these factors, explained much of the variance present in our long-term river water dataset. For example, major patterns of depletion in precipitation isotopes with increasing latitude and elevation were carried through into predictions for rivers and passed downstream using the REC river network (Fig. 6a and d). However, systematic bias was apparent; when compared to measurements from river monitoring sites, the spatiotemporal water balance model showed a tendency to underpredict river water isotope values at the enriched end of the range of validation measurements and overpredict at the depleted end of this range (Fig. 2).

Figure 7 Comparison of modelled δ²H values of river water across two neighbouring North Island basins (Tarawera and Rangitaiki basins). Higher δ²H values from ²H-enriched precipitation lower in these catchments are represented by the water balance model (original). Residual correction using ordinary kriging (OK) extends the enrichment effect of Lake Tarawera (top left of each panel) across the entire region. Using regression kriging (RK), the lake enrichment is mostly confined to downstream reaches. Points denote the locations and observed δ²H values at the gauging sites.

3.6.2 Regression kriging residual corrections improve models of δ²H and δ¹⁸O values in rivers

Across the rivers of New Zealand, a water balance model corrected using the regression kriging of residuals provided improved fits to measurements, compared to a model that used the ordinary kriging interpolation of residuals between validation sites. Patterns of residuals across the validation site network showed distinct patterns, which we attribute to orographic rainfall effects (Kerr et al., 2015; Purdie et al., 2010) and the position of lakes and wetlands in upstream catchments (Yang et al., 2020; Halder et al., 2015). The effect of lakes and wetlands is an example of how landscape processes control δ²H and δ¹⁸O values in rivers; three factors combine to make regression kriging a particularly appropriate method to represent these processes. First, lakes (including artificial lakes and dams) do not cluster predictably across New Zealand, and second, their effects are confined to downstream reaches. Thus, a modelling approach that considers the dendritic nature of river networks in this correction step is likely to better account for this process than one which corrects based only on Euclidean distance (Brennan et al., 2016). Using a term for upstream lake and wetland areas in reach-scale regression kriging incorporates the enrichment effects of these waterbodies and transfers this effect downstream while allowing for mixing. Third, distances between validation sites in the NRWQN network are great, relative to the likely downstream extent of enrichment effects from lakes and wetlands. This raises the likelihood that corrections based only on distance weighting will extend outside of the extent to where their true effect will take place. Representation of evaporation effects in the final ordinary-kriging-corrected and regression-kriging-corrected models is shown using the example of the Lake Tarawera region in Fig. 7. River water observed in the NRWQN validation dataset at the lake outflow is enriched in δ²H and δ¹⁸O relative to upstream precipitation and spatiotemporal model predictions, which do not account for evaporation. Using ordinary kriging, the correction factor based on the residual at the lake outflow site is applied to all nearby reaches and extends outside the catchment boundaries. Using regression kriging correction, the residual correction follows the dendritic nature of the river system and is passed only to downstream reaches.

Both the ordinary kriging and regression kriging approaches to residual correction resulted in a partial correction of the bias associated with orographic rainfall effects; however, spatial extents of the corrections in leeward New Zealand differ between the two kriging methods (Fig. 4). This reflects the data available to predict residuals for both methods; while the ordinary kriging method in limited by the distribution of (and lack of) validation sites across high-elevation areas of southern New Zealand, using regression kriging the correction could be applied at the reach scale via factors that co-vary with orographic effects (e.g. annual rainfall variance and catchment elevation; Table 2). Sharp changes in isotope values of small rivers (Kerr et al., 2015) and snowfall (Purdie et al., 2010) observed across the main divide of the Southern Alps have been associated with orographic effects. That the regression kriging correction method can apply residuals using geographical and climatological gradients that co-vary with orographic weather patterns may partly explain the improved fits of this model to river water measurements nationally.

3.6.3 Regression kriging is vulnerable to overfitting and extrapolation

As mentioned in Sect. 3.3, there are more than five environmental factors in the river network datasets that are spatially correlated with patterns of isotope residuals. However, it is dangerous to include all sensitive factors in regression kriging, as they may cause the problem of overfitting, which is very common in regression analysis and machine learning. Machine learning methods often uses cross-validation techniques to solve the overfitting problem. In this study, there are only 58 sites, which is insufficient for the application of either of data-intensive machine learning methods or nonlinear regression, which increase the number of regression coefficients to fit. Instead, we focused on the five most sensitive environmental factors, following the “one in 10 rule”. Our results indicate that this rule was appropriate for improving model performance and also avoided overfitting; the addition of one more predictor in the regression model caused a wider and unreasonable range of isotope estimates in subsequent national predictions. Extrapolation was another problem we encountered in preliminary regression analysis, as the training dataset did not cover the entire range of values for all catchment characteristics. To nullify this effect, we excluded reaches from the final model where the environmental characteristics of those reaches were outside of the range of the training dataset (e.g. Strahler order 1 streams with upstream coverage of lakes and wetlands greater than 30 % of the total catchment area). Solutions to overfitting and extrapolation could include adding isotope sampling sites and increasing the range of environmental conditions across sampling sites (e.g. setting up river monitoring sites with upstream wetland areas approaching 1) or, for extrapolation, removing reaches from the final isoscape which have values for predictor variables substantially outside the range present across validation sites.

3.6.4 Improved understanding of patterns of precipitation isotope values can improve understanding of hydrological processes across landscapes

Spatial patterns of residuals and predictors of residuals from the spatiotemporal water balance modelling approach could potentially be applied as an investigative tool to describe flow pathways or evaporation processes (Bowen et al., 2011). Regression kriging of these residuals allowed us to identify a suite of catchment environmental variables that were the best predictors of spatial patterns of residuals. However, as described above, patterns of residuals shown in our study are driven to some extent by the systematic bias in high-elevation precipitation isotope input data. Various authors have observed that the origins of air masses contributing to precipitation in New Zealand, and their transit across mountainous regions, contribute strongly to spatial patterns of isotopes in precipitation (McDonnell, 1988; Stewart et al., 1983; Baisden et al., 2016; Kerr et al., 2015). These factors are difficult to capture with empirical models of precipitation isotopes of the kind used in this study (Bowen, 2010); however, potential improvements to the precipitation isotope model could be made by incorporating high-frequency and event-based sampling, by using machine learning prediction methods, or, more simply, by adjusting linear model predictions based on spatial patterns of residuals (Bowen and Revenaugh, 2003; Nelson et al., 2021). Improved regional precipitation isotope input data would raise the visibility of hydrological fluxes (such as high-elevation-derived groundwater contributions to rivers) in the regression kriging correction steps of our method. Similarly, improved regional precipitation isotope models are likely to improve the performance of process-based isotope hydrology models, such as the isotope-enabled coupled catchment–lake water balance model (Belachew et al., 2016) that are designed to quantify hydrological fluxes (e.g. between atmosphere, rivers, lakes, and groundwater) using water isotope data.

4 Final conclusions and implications

We have successfully produced isoscapes of δ²H and δ¹⁸O for the river network of New Zealand by employing a water-balance-based model, enhanced by a regression kriging method to capture catchment characteristics not accounted for by the simple water balance. Although the final model performs well, future research is needed to improve model input data and identify the processes underlying spatial patterns of residuals identified by regression kriging. Further inclusion of these processes which contribute to fractionation and source mixing will reduce the need for correction methods, which are descriptive of current patterns but may not have predictive value with changing climate and hydrological modification.

Distinct differences in spatial patterns of δ²H and δ¹⁸O in river waters to those of precipitation highlight the value of river isoscapes in cross-disciplinary research. Differences between precipitation and river water isotope values were particularly evident at low elevations; New Zealand's high central mountain ranges create strong elevation gradients in precipitation isotopes, and lowland regions receiving more isotopically enriched rainfall are intersected by alpine-fed rivers bearing isotopically depleted water from high elevations. In addition to the hydrological implications described above, quantification of isotopic values for water sources across elevation gradients may be of particular benefit to those studies that attempt to attribute organic material (such as sediment-bound organic material transported in rivers) to particular elevation bands (Upadhayay et al., 2017; Feakins et al., 2016). The river water isoscapes shown in this study are also likely to be appropriate for studies utilising δ¹⁸O values in human tissue to determine historical migration patterns (e.g. King et al., 2021); local drinking water provides a good proxy for the majority of total water intake in humans (Guelinckx et al., 2016; Ehleringer et al., 2008). The development of reach-specific predictions for water stable isotopes may also benefit ecological studies that examine the movement of aquatic organisms (e.g. Brennan et al., 2015; Soto et al., 2013); a fine-scale study of animal movements across river networks can benefit from maps of river water isotope variation across similarly fine spatial scales.

Appendix A

Figure A1 Comparison of measured and modelled water δ²H values across a west–east transect through New Zealand's Southern Alps. Measured values are river samples, taken from small rivers by Kerr et al. (2015), and modelled values are those for precipitation at the nearest virtual climate station network (VCSN) point, based on the coefficients of Baisden et al. (2016), and uncorrected and residual-corrected river model values from this study. All model values show an underprediction of δ²H values at the western (windward) end of the transect and an overprediction at the eastern (leeward) extent of the transect. The RMSEs for the four models across this transect are 13.8 ‰, 12.5 ‰, 9.8 ‰, and 9.6 ‰ for the precipitation model, uncorrected spatiotemporal river model, ordinary-kriging-corrected river model, and regression-kriging-corrected river model, respectively.

Figure A2 Spatial distribution of the five environmental variables used in regression kriging of residuals (Eq. 3). Except for SiteElev, all other factors are averaged across the upstream catchments.

Table A1 Full list of catchment environmental characteristic variables.