Precipitation isoscapes for New Zealand: enhanced temporal detail using precipitation-weighted daily climatology†

ABSTRACT Predictive understanding of precipitation δ2H and δ18O in New Zealand faces unique challenges, including high spatial variability in precipitation amounts, alternation between subtropical and sub-Antarctic precipitation sources, and a compressed latitudinal range of 34 to 47 °S. To map the precipitation isotope ratios across New Zealand, three years of integrated monthly precipitation samples were acquired from >50 stations. Conventional mean-annual precipitation δ2H and δ18O maps were produced by regressions using geographic and annual climate variables. Incomplete data and short-term variation in climate and precipitation sources limited the utility of this approach. We overcome these difficulties by calculating precipitation-weighted monthly climate parameters using national 5-km-gridded daily climate data. This data plus geographic variables were regressed to predict δ2H, δ18O, and d-excess at all sites. The procedure yields statistically-valid predictions of the isotope composition of precipitation (long-term average root mean square error (RMSE) for δ18O = 0.6 ‰; δ2H = 5.5 ‰); and monthly RMSE δ18O = 1.9 ‰, δ2H = 16 ‰. This approach has substantial benefits for studies that require the isotope composition of precipitation during specific time intervals, and may be further improved by comparison to daily and event-based precipitation samples as well as the use of back-trajectory calculations.


Introduction
Isoscapes, or predictive modelling of isotope ratios across landscapes from local to global scales, have become an important paradigm in the application of isotope science to environmental problems. Precipitation isoscapes represent the δ 2 H and δ 18 O of precipitation, and therefore provide essential inputs to any hydrological or biological or environmental application of H and O isotopes [1]. The successful methodologies reported by Bowen et al. [1] for calculating precipitation isoscapes rely on mean annual or mean monthly precipitation amount weighted δ 2 H and δ 18 O. Below, we outline a wide variety of applications utilizing spatial and temporal information represented by precipitation isoscapes. For many examples, we identify that improved temporal resolution would be beneficial in locations such as New Zealand, in order to capture weekly to seasonal variation in climate due to blocking [2,3], synoptic properties of rain-inducing air masses [4], with an ultimate goal of understanding sub-hourly variation during precipitation events [5].
Understanding the potential need for enhanced temporal resolution in isoscapes must be placed within the context of various applications that will benefit from the improved understanding of New Zealand's variable exposure to alternating sub-Antarctic and subtropical air masses. Historic applications in New Zealand have centred on hydrology, where hydrograph separation and study of runoff generation processes during storms benefited from the large and predictable variation in water isotopes between storms of varying moisture source [6]. Due to the importance of variation in δ 2 H and δ 18 O during storm events, these applications are among the most demanding of high temporal resolution [5]. Building on the foundations of isotope hydrology, strong potential exists for δ 2 H and/or δ 18 O to provide highly detailed paleoclimate information. Currently, isotopic records of paleoclimate for New Zealand consist mainly of δ 18 O in speleothems, which show a somewhat ambiguous relationship with temperature [7]. New areas with potential for much higher resolution records include tree-ring records [8], which can potentially be interpreted at a sub-annual resolution [9]. These proxies may be augmented by compound-specific δ 2 H in plant waxes or δ 18 O of cellulose in sediments and peats, respectively.
The physiological understanding of H and O isotope systematics from precipitation, to water sources or soil, and ultimately to living organisms also underpins a range of other applications. Forensic science broadly represents one range of these applications [1,9]. Forensic methodologies have now been applied to materials such as feathers [10], which can define the provenance of indigenous artifacts, or delineate the migration patterns of birds. Similarly, it has been shown that the source location of New Zealand's leading export, whole milk powder, can be determined using bulk or compound-specific δ 2 H assays [11]. Finally, New Zealand's geographic isolation has defined biosecurity as a key application, where interest has focused on the ability to identify whether a single captured insect is derived from a locally established population, or of overseas origin [12].
In each of the physiological applications described above, integrated δ 2 H and δ 18 O is required for a specific time intervalwith duration of weeks or monthsduring which tissues have formed. To address problems such as these, we briefly describe a dataset of δ 2 H and δ 18 O from monthly meteoric water samples collected for 3 years at over 50 sites throughout New Zealand. We then develop a conventional mean annual isoscape model as well as a model using weighted daily climate variables designed to provide a new isoscape with spatial and temporal resolution appropriate for biological applications for New Zealand.

Measurements
A dataset of the spatial distribution of the H and O isotopic composition of New Zealand rainfall was collected in a separate study [13]. Briefly, 51 rainwater collection sites were selected to sample the range of environments found in New Zealand. Monthly sampling was conducted from January 2007 to December 2009. Most sites were operated by volunteers, necessitating that some sites were dropped over time and replaced to maintain similar geographical coverage. Samples were collated and measured at the Iso-trace Laboratory, University of Otago. For δ 2 H, small (0.3 µL) aliquots of water were reduced to H 2 gas over hot glassy carbon using a temperature conversion elemental analyzer. The isotopic composition of the resulting gas was measured via continuous flow on a Delta V mass spectrometer (Thermo Finnigan, Bremen, Germany). Similarly, δ 18 O was measured using a Thermo Finnigan GasBench and Delta Plus Advantage continuous flow mass spectrometer. Analyses were performed on helium containing 0.3 % carbon dioxide that had equilibrated with 0.5 mL of the water sample at a constant temperature of 25°C for 18 h. Data were calibrated against three internal laboratory standards (SEA TAP, and ICE, spanning a range of δ 2 H, from −1.5 to -262.7 ‰) and reported with respect to VSMOW-SLAP scales.
A subset of samples was measured for quality control at the National Isotope Centre, GNS Science, Lower Hutt, New Zealand, using the following methods: The δ 2 H values of water were obtained by reduction on a GVI PyrOH attached to a GVI IsoPrime mass spectrometer using 5 µL of water injected into a helium stream through a quartz reactor filled with chromium granules and quartz chips at 1050°C. The resulting H 2 was analysed in continuous flow mode on the IsoPrime. Similarly, δ 18 O was obtained through CO 2 equilibration in Exetainer™ vials on a Gilson 222XL, and measurement via a dual-inlet IsoPrime mass spectrometer. Both isotopes were reported with respect to VSMOW-SLAP scales, normalized to internal standards: INS11, INS9, and MM1. Subsequent validation data were also analysed at GNS Science, using the mass spectrometers described above or Los Gatos Research LWIA and IWA-35 laser spectrometers, to INS11 and CM1 internal standards.

Model
Monthly precipitation from the network of 51 sites across New Zealand was measured for δ 2 H and δ 18 O during the period 2007-2009 [13]. The measurements were used to develop a multiple linear regression model to predict daily δ 2 H and δ 18 O values across New Zealand from two geographic variables and five daily climate variables from the Virtual Climate Station Network (VCSN) using JMP 8 software (SAS Institute, Cary, NC, USA). The VCSN is a climate database, developed and maintained by the National Institute of Water and Atmospheric Research (NIWA), that uses interpolation techniques to create a daily climate record at every point on a ∼5-km grid covering all of New Zealand [14]. To maximize the representativeness of output statistics, the closest VCSN data were used for all precipitation collection sites, regardless of other potential climate data sources. The mean annual δ 2 H isoscape reported here was compiled using ordinary least squares regression in SPSS 19.0 (SPSS, Chicago, IL, USA). A total of 51 sites with at least one complete year of data were used in the regression, 47 of which had ≥2 years of data. To provide an isoscape representing a best estimate of average climate, results were mapped using average climate data for the years 1950-2000, mapped at 2.5 arcmin resolution, and available from http://www.worldclim.org as reported in Rogers et al. [10].
To obtain an isoscape that could be considered statistically valid on a daily timescale, daily climate variables were weighted by daily precipitation amounts to generate monthly weighted climate means corresponding to precipitation samples measured for δ 2 H and δ 18 O. The rationale for the approach is that weighting places emphasis according to climate parameters when rainfall occurs, and not average climate for the month. Thus, results of the regression will reflect the few days each month when substantial rainfall occurred. Technical details of the procedure are as follows: Precipitation itself was not used as an independent predictor to avoid expected instability in model predictions when precipitation events of the same size and duration were split over two days, versus one single day. The regression was compiled in JMP 8 software (SAS Institute) using stepwise backward regression, supervised to eliminate highly correlated variables and use a nearly identical independent variable set for both δ 2 H and δ 18 O. Backwards elimination was justified because the significance of many parameters decreased markedly when other parameters were excluded from the regression. Plots showing the 'leverage' of each variable were used to check for and remove parameters with high levels of colinearity.
The weighted means of climate variables were regressed as independent variables explaining all available monthly δ 2 H and δ 18 O. Deuterium excess, d, was also regressed as a dependent variable. The mean annual δ 2 H, δ 18 O, and d were evaluated for sites with 2 or more years of data by regressing the precipitation-weighted average of monthly predictions versus measurements. For comparison to these station δ 2 H averages, the model was applied to data from 2007 to 2009 to generate map-based estimates of the δ 2 H values of precipitation, in Matlab (Natick, MA, USA). Table 1 shows that a regression using elevation, temperature, and precipitation yields a significant and predictive regression of mean annual δ 2 H values based on 51 stations across New Zealand. The explanatory power of the model is high (R 2 = 0.86), compared to similar regressions presented by Bowen et al. [1] for global, United States, and Mediterranean region datasets (R 2 in range 0.58-0.79). No additional benefits were obtained by including a term in the regression to account for spatial autocorrelation, confirming that climate and elevation are the best explanatory variables available. Additional details and a map derived from the parameters in Table 1 appear in Rogers et al. [10].

Results and discussion
Availability of a series of daily VCSN variables for New Zealand enables more temporal details to be considered and evaluated as a tool for isoscape development. The advantage of daily climate data is that weighting each climate parameter by daily rainfall allows the regression to focus appropriately on climate when the rainfall occurred, for each monthly sample. This procedure requires two steps to enable comparison to a mean annual isoscape, described using a model such as the regression reported in Table 1. First, the monthly precipitation-weighted average of the daily climate variables is regressed against the monthly precipitation δ 2 H and δ 18 O values (Table 2). Then, the precipitation-weighted average δ 2 H and δ 18 O values from the monthly model and the monthly measurements at each station can be assessed using regression (Figure 1). The advantage of this procedure is that the full dimensionality of the original data is preserved in primary (monthly) regression. The larger dimensionality of the data allowed a more complex set of climate variables to be included in the regression, using a backwards elimination procedure. In addition, all stations can be included without introducing bias, regardless of missing samples and incomplete data.
For monthly δ 2 H predictions, the model yielded R 2 = 0.43 and RMSE = 16 ‰ ( Table 2, Supplemental Figure S1(a)), and R 2 = 0.79 and RMSE = 5.5 ‰ for long-term average predictions (Figure 1(a)). For δ 18 O, the model yielded R 2 = 0.46 and RMSE = 1.9 ‰ for monthly δ 18 O predictions (Table 2, Supplemental Figure S1(b)) and R 2 = 0.85 and RMSE = 0.6 ‰ for long-term average predictions (Figure 1(b)). Due to the scarcity of published predictive models for d-excess, it is also useful to report the development of a predictive model for d-excess using this methodology. For monthly predictions, the model yielded only R 2 = 0.21 and RMSE = 4.0 ‰ ( Table 3, Supplemental Figure S1(c)). However, long-term average predictions yielded R 2 = 0.89 and RMSE = 0.9 ‰ (Figure 1(c)). The RMSE values improved from individual monthly predictions, to long-term average predictions, by factors of 2.9, 3.1, and 4.4 for δ 2 H, δ 18 O, and d, respectively, indicating that the measurements are largely independent of one another, resulting in uncertainty reduction during averaging.
The predictive models for δ 2 H and δ 18 O contained similar parameters, but with some differences. Elevation, latitude, mean sea level pressure, solar radiation, minimum temperature (T min ), and wind speed had consistent signs, and magnitudes with coefficients for δ 2 H approximately an order of magnitude greater than those for δ 18 O, as expected. For δ 2 H, the elevation coefficient (-0.017 ± 0.002) was similar to the mean annual model in Table 1 (-0.013 ± 0.003). The latitude coefficients were similar to those in Bowen et al. [1]. The δ 2 H regression contained soil temperature at 10 cm, which was excluded from the δ 18 O model due to the lack of significance. Similarly, the δ 18 O regression found two additional terms, potential evapotranspiration (PET) and vapour pressure, to be highly significant, but these terms were not included in the δ 2 H model. Due to the variation in d within the dataset, it was not surprising that the terms varied in significance between δ 2 H and δ 18 O regressions, especially vapour pressure which has long been known as a control on kinetic fractionation explaining d [15].
In this regard, it was unsurprising that the predictive model for d contained vapour pressure and relative humidity, as well as latitude and soil temperature at 10 cm. The inclusion of both vapour pressure and relative humidity can be presumed to relate to the fact that these related variables differ in their formulation and in the degree to which they include temperature, and may have linear relationships with d for different reasons. Terms such as T min and PET, which also differed in sign between the δ 2 H and δ 18 O regressions, and soil temperature at 10 cm in the regression for d, may relate to evaporation (or transpiration) fluxes from the land surface, as causes of variation in d. Overall, the different terms in the δ 2 H, δ 18 O, and d regressions appear reasonable, and if tested in other regions, and formulations with increasing mechanistic understanding, could add nuances to our understanding of isotope systematics in precipitation.
A δ 18 O isoscape resulting from our model for 2007-2009 is shown in Figure 2, with comparison to precipitation-weighted average measurements (black circles) for each station. This visualisation augments Figure 1(c) with an alternate representation of the model's fidelity. Most stations are in excellent agreement with the surrounding isoscape, and those that differ are largely in areas of steep mountainous topographyareas of challenges and opportunities for improving high-resolution isoscapes [16,17]. Monthly data overlain onto the modelled isoscape can also be viewed as a video file (http://j.mp/ NZisoprecip-mo), revealing greater model-data disagreement on a monthly time step. Daily isoscapes can also be visualised (http://j.mp/NZisoprecip-da), noting that measurements are not available and cannot be overlain onto the model. The isoscape model can also be applied to long-term climate data, creating a long-term average more representative than is possible using 2-3 years of isotope data.
The approach has been validated through an application of the δ 2 H isoscape: modelled precipitation δ 2 H values during August-December 2009 were strongly related to the bulk δ 2 H value of whole milk powder (R 2 = 0.86, RMSE = 3 ‰), as well as the compound-specific δ 2 H of C14:0, C16:0, C18:1 fatty acids [11]. In the case of transfer of the δ 2 H signal from precipitation, to soil water and then pasture grass, and finally milk, a 5-month average appeared to be a sufficient averaging period for the temporally explicit isoscape method to yield very useful results.
A greater challenge is the transfer of precipitation into stream water discharge, where it is desirable to accurately estimate both the amount and isotopic composition associated with each precipitation event. In a region where no precipitation collection was undertaken, event and monthly bulk precipitation samples were collected near Carterton, New Zealand, from October 2010 until January 2013. Stream water was also collected from three locations in the Mangatarere River.
The results (Supplemental Figure S2) illustrate that the monthly collections were often well represented by the isoscape prediction model, but individual precipitation events were poorly represented in both size and isotope ratio. The results support the view that the temporally resolved model produces daily results that are statistically valid, but with a high level of variability limiting immediate use of model predictions at daily resolution. The visualization in Supplemental Figure S2 provides an understanding of the seasonal and event-driven amplitude of variation in precipitation δ 2 H, and its relationship to the much smaller variation in streamwater. The large variation suggests that the following additional steps to improve predictions at the event scale could be beneficial. Our model methodology encourages the collection of daily or event samples, because these can be included in the regression model alongside monthly data, as long as the weighting procedures are appropriate. There is reason to expect that further progress will be enabled by daily or event collections at a limited number of stations, where field collection and analytical costs can be optimised within acceptable bounds. In addition to allowing daily or event-based data to be utilised alongside monthly data, and reducing issues associated with missing data, it will become practical to use new independent data such as air mass back-trajectory information [18]. Higher temporal resolution data from within events will also be valuable, better enabling inclusion of precipitation amount in the predictive regression, and enabling progress on precipitation from air masses crossing mountain terrain.

Conclusions
A method for generating precipitation δ 2 H and δ 18 O isoscapes using precipitation-weighted daily climate yields similar mean annual isoscapes to conventional methods. Advantages of the proposed method are the ability to utilise the full dimensionality of the dataset, including stations with missing data, and the ability to make statistically valid predictions for specific periods based on actual climate data. The new method allows visualization of isoscapes at monthly or daily time steps, but achieves most useful accuracy when integrating across several months. At longer timescales, the new method also appears to produce accurate integrated predictions of d-excess. Predictions have been validated in a study of bulk and compound-specific δ 2 H in whole milk powder and associated fatty acids. The method opens room for additional progress through daily, event-based, or intra-event sampling of precipitation isotopes, inclusion of back-trajectory information or modelling, and better inclusion of precipitation amounts. This approach has substantial benefits for studies such as catchment hydrology, tracing animal migration, agricultural authentication, biosecurity, and forensic studies that may require the isotope composition of precipitation during specific periods.