Response: We have interpreted this comment as a summary of the bullet points given below under the reviewer’s section (a), and address these individually.

R: We agree that our stream measurements dataset, and modelled reach scale values (i.e. those available at https://shiny.niwa.co.nz/nzrivermaps/ as shown below)  may be useful contributions regionally.

However, we note that the method used to derive modelled reach scale values improves upon a method originally applied to the continental USA (Bowen et al. 2011). Our improved method is transferrable to other parts of the world and may thus be of value for an international audience. To better clarify to readers of the novelty of this method and the improvement it offers relative to other methods of mapping river water isotope values, we will improve description of the regression kriging method, its uses elsewhere, and the novelty in its application to river water isoscapes (as suggested by reviewer 2).

R: We think a brief review of the history of environmental isotope studies in New Zealand is a good idea. We will add text to the introduction to supplement the description of the literature validation dataset in the methods (section 2.4). In response also to comments below from both reviewers, will add a table to show the sources, size and duration of datasets and how these data are used in developing, calibrating and validating the models.

Response: We will include the GeoTIFF files for precipitation and river isoscapes in our supplementary materials.

Response: Figure 1 already shows the river network (panel F) and the locations of validation sites (panel G). It is an issue that at national scale, in panel F the smaller reaches and catchments blend together so that only larger rivers are visible. For this reason, panels A-F show only the Canterbury region of New Zealand. We cannot show every catchment in high detail, but a ‘zoomed-in’ example comparing model results for individual reaches to values from monitoring sites is presented in Figure 7.

For the reader with further interest in the river network and catchments we will add more detail as follows:

• Add reference Smith and McBride (1990) to panel G in Figure 1; this reference describes the design of New Zealand’s national river water quality network, and the catchments from which monthly isotope samples were taken.
• Add text at the start of section 3 to give the reader better access to information about monitoring sites and their catchments, including design of the monitoring network (Smith and McBride 1990), and descriptions of physical (catchment), flow and chemical conditions at monitoring sites (Davies-Colley et al. 2011; Julian et al. 2017; Yang et al. 2020).

Response:

Regarding the presentation of results related to surface flow mixing or patterns:

As noted in L. 16, We used a water balance-based method to generate the river isoscape. Patterns of surface flow and mixing are therefore represented by the isoscape outputs in figures 6 and 7, in the reach scale values available at https://shiny.niwa.co.nz/nzrivermaps/   and will be available in the GeoTIFF files for river isoscapes we will add as supplementary material.

Regarding the comment about the main contribution of the work is not being about isotopes in animals or plants:

We will remove mention of animals and plants from the first sentence of the manuscript. This will read: ‘Stable isotope ratio measurements (isotope values) of surface water reflect hydrological pathways, mixing processes, and atmospheric exchange within catchments.’.

We will slightly alter the last sentence of the abstract to make it clear that we haven’t measured any animals or plants. This will read: ‘The resulting river water isoscapes have potential applications in ecological, hydrological and provenance studies for which understanding of spatial variation in surface water isotope values is required' 

Response: We have interpreted these comments as a summary of the bullet points given below under this review section (b), and address these individually.

 Response: We agree. We will address this by providing more background on the differences between ordinary kriging and regression kriging (as requested by reviewer 2) and using the term ‘simple distance-based’ to describe ordinary kriging in lines 15 and 54.

Response: We will provide more detail and support our choice of method with a reference at this point. E.g. ‘From the list of independent variables in Table S1, five were selected for the regression analysis based on BIC (Baysian Information Criteria), following the “one in ten rule” (e.g. Harrell Jr (2015)), i.e. one predictive variable can be included for every ten sites in the dataset.’

We will add t values and P values to this table.

Response: Really good idea. Thanks. Reviewer 2 also had trouble working out which datasets were used in which step, and what data these contained. We will add a table to make this clearer.

Response: We will try to make this clearer using the addition of a table as described in the previous comment.

As described in a response to reviewer 2, the 2017 to 2020 river water monitoring data were monthly samples over three years (36 samples per site) from 58 sites spread across the major catchments of New Zealand. Site mean values of δ²H and δ¹⁸O from this dataset are appropriate for correcting the river isotope model; the model gives an estimate of ‘average’ river water isotope values for each reach in the river network.

Other data collated from the literature used in final checking of the model (figure 5, and supplementary material S1) is largely from ‘one-off’ samples from river reaches, which are less appropriate for correcting the model because river water isotope values vary seasonally, and under changing flow conditions (Yang et al. 2020). These data do however provide a good independent check of how the model performs compared to other model approaches.

Response:

We will make sure this field precipitation dataset and its use in model checking are clearly described using the additional table suggested above.

The dataset of field precipitation samples used for this correlation/model checking contained monthly values from 51 sites between 2007 and 2010. So, ca. 1400 data points for computing the monthly correlation, and 51 data points for the annual average correlation (not 4). We will add these samples size values to the manuscript text alongside the R² and RMSE values. We note that this correlation method and the field dataset are the same as used by Baisden et al. (2016). Using the same corelation method and field dataset allowed us to check our precipitation model replicated their published one well.

Response: We didn’t apply nonlinear regression simply because it would increase the complexity for parameter estimates. This was inappropriate given limited number of sites (58) available for validation. We will make this clear in the manuscript.

Response: We will revise this section for clarity.

Spatial patterns of residuals in our method, and predictors (e.g. those in Table 1) could be used to increase understanding of hydrological processes. A simple example of this is that upstream wetland and lake area, which leads to higher evaporative fractionation and thus higher river water δ²H and δ¹⁸O values, explained spatial patterns of residuals in our study, which used the (Bowen et al. 2011) water balance model that assumes no evaporative fractionation.

A similar approach may be taken for lowland reaches gaining a large portion of their flow from high-elevation-derived groundwater. These reaches may show up as more isotopically negative than would be expected based on recharge and surface routing of local, low-elevation rainfall.

HOWEVER – our manuscript shows that this type of approach, and similar approaches in isotope enabled hydrological models (e.g.  Belachew et al. (2016)) are reliant on the accuracy of the precipitation isotope model. We believe that some of the variables in table 1 reflect correction of spatial inaccuracies of the precipitation model.

We will adjust the discussion accordingly.

Response: Our method is fairly robust in this respect. Because both input data and data used to correct the model are averaged to give ‘steady state’ values, seasonal variation in contributions of surface and groundwater flows (which may have differing isotope values) to rivers is incorporated.

We will add brief discussion on this point.

We have interpreted this comment as a summary of the specific bullet points given below under this review section (c), and address these individually.

We agree that the discussion of the effects of origin of air masses on precipitation isotope values currently looks speculative because we haven’t referenced previous work well enough. We will add references to studies of isotopes in precipitation in New Zealand (e.g. (McDonnell 1988) to back up this point.

We feel our discussion of regional orographic effects is well supported by the work of Purdie et al. (2010) and Kerr et al. (2015), which we have referenced in the manuscript. We have gone to further effort to back up our statements using Appendix Figure 1 and its accompanying text. However, we will add international references showing the same effects in other mountainous regions worldwide to support our statements.

Response: Firstly, we will directly reference Figure 2 in the statement on L. 225-226; i.e. ‘of the eight sites where predicted δ¹⁸O values exceed average measured δ¹⁸O values by > 1‰ (Figure 2), seven are in alpine-fed rivers on the leeward east of New Zealand.

We agree that there is some unnecessary repetition between sections 3.2 and 3.4. We will work to reduce or remove this.

We will revise the discussion to improve the description of links between:

1. The predictors of residuals in Table 1
2. The spatial inaccuracy of the precipitation model (and its likely causes)
3. Our ability to improve understanding of processes in hydrology using our approach.

Put simply, Table 1 currently contains predictors that correct for spatial inaccuracy of the precipitation model. If we can improve the precipitation model, our method will be more useful for understanding of processes in hydrology (such as evaporation and groundwater contributions to surface water) that change in isotope values of river water.

Response: As above, we will revise the discussion to improve the description of links between:

1. The predictors of residuals in Table 1
2. The spatial inaccuracy of the precipitation model (and its likely causes)
3. Our ability to improve understanding of processes in hydrology using our approach.

 Response: Good point. In fact, the ‘usAnRainVar’ variable in Table 1 is strongly correlated with aspect. We will make this clearer as described above.

Response: The reviewer is right that the current maps have potential use in hydrological studies. With the reviewer’s comments and the readership of HESS in mind, we will make modifications to the abstract, introduction and discussion to lessen the focus on ecological implications of this work and increase focus on hydrological uses and implications.

We do not feel that the absence of nitrogen data from our paper negates the usefulness of our work to (for example) ecological research. While having MORE tracers is almost always better in mixing models (Fry 2006), hydrogen and oxygen stable isotopes are useful nonetheless for aquatic ecology (Soto et al. 2013).

We do not agree that the results of this work are unlikely to be useful for the ecological purposes we have outlined in our manuscript. The geographical distributions of hydrogen and oxygen isotopes in precipitation and surface water form underpin a rich and growing body of research into animal migrations, as well as other cross-disciplinary uses. Quoting from Bowen et al. (2009) ‘Isoscapes have great power as a cross-disciplinary research tool, as exemplified by the translation of hydrology-focused GNIP [Global Network of Isotopes in Precipitation] data into tools for animal migration research.’. Examples of ecological (migration) research based on GNIP hydrogen and oxygen isotope data are included in a review by Hobson and Wassenaar (2018). The Global Network for Isotopes in Rivers (GNIR) has similar aims. Quoting from Halder et al. (2015) ‘The aim of the GNIR programme is to collect and disseminate time-series and synoptic collections of riverine isotope data from the world’s rivers and to inform a range of scientific disciplines including hydrology, meteorology and climatology, oceanography, limnology, and aquatic ecology.’ However, the reviewer’s comments make it plain that we have not conveyed this potential for cross-disciplinary use of our work adequately. To address this, we will add brief but specific examples to the manuscript on this topic to section 4.  

 R: We will extend the focus of this paragraph outside of the scope of the current study, towards more general discussion of using river isotope models to understand hydrological processes.

We will restructure this paragraph as follows:

1. Some isotope-enabled hydrological models (e.g. Belachew et al. (2016)) use precipitation isotope models as input data to give improved estimates of fluxes between components of the hydrological cycle.
2. The accuracy of these flux estimates relies partly on accuracy in input data from precipitation isotope models
3. Data from precipitation isotope models will always be imperfect, but improvements in the accuracy of precipitation isotope models can improve our understanding of flow pathways and evaporation processes at landscape scales.

Response:

We will add t values and p values to Table 1 to better support the discussion around upstream lake and wetland area effects on δ²H and δ¹⁸O values of river water. We will refer the reader to these results in section 3.3.

Open water behind dams is included in the variable usLWArea. We will add a specific reference to dams to section 3.6.

For clarity, we will replace abbreviated variable names (e.g. ‘usLWArea’) in the results and discussion text with full variable names (e.g. ‘Upstream lake and wetland area’).

Line 16 states ‘We used a water balance-based method, which represents patterns of surface flow and mixing’. Thus, mixing results are incorporated into the water balance results shown in (for example) Figure 2, 6 and 7, and in the online maps shown below.  

Response:

We will add t values and p values to Table 1 to better support this discussion section.

Discussion of relationships between environmental variables and river water δ²H and δ¹⁸O in Ln355 - Ln379 is backed up by multivariable regression results. Importance ranks for this regression for δ²H and δ¹⁸O residuals are already given in table 1, based on t statistics. This regression used δ²H and δ¹⁸O data from across New Zealand, not just the catchment in figure 7. Our discussion on Ln355 - Ln379 is limited to statistically significant predictors shown in table 1, of which upsteam lake and wetland area is one.

We could produce plots similar to Figure 7 for all the main catchments in New Zealand, but it is not practical to show them all in the manuscript. Figure 7 gives an example. We have provided access to data shown in Figure 7 at https://shiny.niwa.co.nz/nzrivermaps/. A South Island example is shown below:

Monitoring data (i.e. equivalent to the points in Figure 7, but across major catchments nationally) are available via the IAEA WISER portal.

Response: Certainly. We will add some detail on this to the section in L. 364-378.

 Response: Really good point. We will add more discussion to section 4, focussing on implications of differences in isotopes in precipitation and those shown in our runoff maps. In terms of quantification, to some degree this is already visible in the isotopic differences between rivers fed by high elevation recharge and those fed by local lowland recharge (see above). We will add some text to this effect.

Thank you. As above, to better clarify to readers the international transferability of our work, we will improve description of the regression kriging method, its uses elsewhere, and the novelty in its application to mapping river water isotope values (as also suggested by reviewer 2).

References:

Baisden, W.T., E.D. Keller, R. Van Hale, R.D. Frew, and L.I. Wassenaar. 2016. Precipitation isoscapes for New Zealand: enhanced temporal detail using precipitation-weighted daily climatology. Isotopes in Environmental and Health Studies 52: 343-352.

Belachew, D.L., G. Leavesley, O. David, D. Patterson, P. Aggarwal, L. Araguas, S. Terzer, and J. Carlson. 2016. IAEA Isotope-enabled coupled catchment–lake water balance model, IWBMIso: description and validation. Isotopes in Environmental and Health Studies 52: 427-442.

Bowen, G.J., C.D. Kennedy, Z. Liu, and J. Stalker. 2011. Water balance model for mean annual hydrogen and oxygen isotope distributions in surface waters of the contiguous United States. Journal of Geophysical Research: Biogeosciences 116.

Bowen, G.J., J.B. West, B.H. Vaughn, T.E. Dawson, J.R. Ehleringer, M.L. Fogel, K. Hobson, J. Hoogewerff, C. Kendall, and C.T. Lai. 2009. Isoscapes to address large‐scale earth science challenges. EOS, Transactions American Geophysical Union 90: 109-110.

Davies-Colley, R.J., D.G. Smith, R.C. Ward, G.G. Bryers, G.B. McBride, J.M. Quinn, and M.R. Scarsbrook. 2011. Twenty Years of New Zealand’s National Rivers Water Quality Network: Benefits of Careful Design and Consistent Operation1. JAWRA Journal of the American Water Resources Association 47: 750-771.

Fry, B. 2006. Stable isotope ecology: Springer.

Halder, J., S. Terzer, L. Wassenaar, L. Araguás-Araguás, and P. Aggarwal. 2015. The Global Network of Isotopes in Rivers (GNIR): integration of water isotopes in watershed observation and riverine research. Hydrology and Earth System Sciences 19: 3419-3431.

Hobson, K.A., and L.I. Wassenaar. 2018. Tracking animal migration with stable isotopes: Academic Press.

Julian, J.P., K.M. de Beurs, B. Owsley, R.J. Davies-Colley, and A.G.E. Ausseil. 2017. River water quality changes in New Zealand over 26 years: response to land use intensity. Hydrol. Earth Syst. Sci. 21: 1149-1171.

Kerr, T., M. Srinivasan, and J. Rutherford. 2015. Stable water isotopes across a transect of the Southern Alps, New Zealand. Journal of Hydrometeorology 16: 702-715.

McDonnell, J.J. 1988. The age, origin and pathway of subsurface stormflow in a steep humid headwater catchment. PhD thesis, University of Canterbury Canterbury, New Zealand.

Purdie, H., N. Bertler, A. Mackintosh, J. Baker, and R. Rhodes. 2010. Isotopic and elemental changes in winter snow accumulation on glaciers in the Southern Alps of New Zealand. Journal of climate 23: 4737-4749.

Smith, D.G., and G.B. McBride. 1990. New Zealand's national water quality monitoring network - design and first year's operation. JAWRA Journal of the American Water Resources Association 26: 767-775.

Soto, D.X., L.I. Wassenaar, and K.A. Hobson. 2013. Stable hydrogen and oxygen isotopes in aquatic food webs are tracers of diet and provenance. Functional Ecology 27: 535-543.

Yang, J., B.D. Dudley, K. Montgomery, and W. Hodgetts. 2020. Characterizing spatial and temporal variation in 18O and 2H content of New Zealand river water for better understanding of hydrologic processes. Hydrological Processes.

Responses to referee comment R2 on:

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

By Bruce D. Dudley, Jing Yang, Ude Shankar and Scott Graham

 

Response: We agree that the majority of the variation in river water isotope values across New Zealand can be explained by the water balance model used by Bowen et al. (2011) for the continental USA. This is a good point, and we will change the abstract and conclusions to make it clear. 

Nevertheless, the additional effort put into regression kriging does result in improved predictions, and there are some areas of river networks, such as downstream of dams and wetlands, where it appears particularly beneficial.

Response: Yes, regression kriging is tool used elsewhere, particularly in soil mapping (e.g. see (Hengl et al., 2007; Keskin and Grunwald, 2018). We will add text to the introduction to describe other common uses of regression kriging and discuss the novelty (and appropriateness) of its application to this particular problem. 

Response: Yes, while the uncorrected river model was generated using 20 years of gridded climate data, only a 3-year period of river isotope data were available to correct the river model. We will add text (and a table) to the manuscript to make this clear. 

Because the uncorrected river model gives an estimate of flow-weighted average δ²H and δ¹⁸O at any reach, we needed comparable average values at each of our 58 river monitoring sites to correct it. While more data is always better, annual means (of monthly) δ²H and δ¹⁸O values at our monitoring sites were relatively consistent across the three years of monitoring data we have, so we are confident that our 2017-2020 dataset is adequate for the purposes of model correction. As a rough illustration of this consistency between years, the isotope biplots below show all monthly river water samples from our 58 validation sites, with additional years of data showing in red, then blue. We note also that the monitoring dataset we have collected is the largest and by far the longest record of river water isotopes available for New Zealand.

In response to the reviewer’s question about expectations for the future: we will make our code publicly accessible so that updates to modelling methods, and additional data can iteratively improve these maps. We continue river sampling at monitoring sites, and we are also now working to improve the national precipitation isotope model. These new input and validation data will be incorporated into revised maps. When available, updated maps will be made available online via https://shiny.niwa.co.nz/nzrivermaps/ 

Additional (measured) river data will be added to the dataset we have provided to the IAEA WISER database.  

Response: This period of analysis was again limited by availability of field (measured) data. We will add text to the manuscript to make this clear, and, as suggested by reviewer 1, a table showing the different datasets used in model development, correction and checking. The reviewer has also highlighted a typo – the dataset of Baisden et al. (2016) extended into early 2010 at some sites. We will correct this. 

References:

Bowen, G. J., Kennedy, C. D., Liu, Z., and Stalker, J.: Water balance model for mean annual hydrogen and oxygen isotope distributions in surface waters of the contiguous United States, Journal of Geophysical Research: Biogeosciences, 116, 2011.

Hengl, T., Heuvelink, G. B. M., and Rossiter, D. G.: About regression-kriging: From equations to case studies, Computers & Geosciences, 33, 1301-1315, https://doi.org/10.1016/j.cageo.2007.05.001, 2007.

Keskin, H. and Grunwald, S.: Regression kriging as a workhorse in the digital soil mapper's toolbox, Geoderma, 326, 22-41, https://doi.org/10.1016/j.geoderma.2018.04.004, 2018.