Global isoscapes for δ 18 O and δ 2 H in precipitation : improved prediction using regionalized climatic regression models

A regionalized cluster-based water isotope prediction (RCWIP) approach, based on the Global Network of Isotopes in Precipitation (GNIP), was demonstrated for the purposes of predicting pointand large-scale spatio-temporal patterns of the stable isotope composition ( δ2H, δ18O) of precipitation around the world. Unlike earlier global domain and fixed regressor models, RCWIP predefined 36 climatic cluster domains and tested all model combinations from an array of climatic and spatial regressor variables to obtain the best predictive approach to each cluster domain, as indicated by root-mean-squared error (RMSE) and variogram analysis. Fuzzy membership fractions were thereafter used as the weights to seamlessly amalgamate results of the optimized climatic zone prediction models into a single predictive mapping product, such as global or regional amount-weighted mean annual, mean monthly, or growingseasonδ18O/δ2H in precipitation. Comparative tests revealed the RCWIP approach outperformed classical global-fixed regression–interpolation-based models more than 67 % of the time, and clearly improved upon predictive accuracy and precision. All RCWIP isotope mapping products are available as gridded GeoTIFF files from the IAEA website (www.iaea.org/water ) and are for use in hydrology, climatology, food authenticity, ecology, and forensics.


Introduction
Spatial patterns in stable hydrogen and oxygen isotope ratios of precipitation were first observed in the 1950s (Epstein and Mayeda, 1953;Friedman, 1953;Dansgaard, 1954;Craig, 1961), and increasingly revealed as long-term δ 2 H and δ 18 O data sets from around the world accumulated in the Interna-tional Atomic Energy Agency's (IAEA) global network of isotopes in precipitation (GNIP -Fig. 1) (Dansgaard, 1964;Rozanski et al., 1993;Aggarwal et al., 2010;IAEA/WMO, 2013).The strong covariance of δ 2 H and δ 18 O in precipitation (Craig, 1961) supported the construction of local and regional isotopic "meteoric water lines", which provided the basis, for example, upon which to assess the origin of modern and ancient ground water (Rozanski, 1985;Clark and Fritz, 1997) and its interaction with surface water resources (Kendall and McDonnell, 1998;Froehlich et al., 2005).
The past decade has experienced increasing demand for accurate spatio-temporal predictions of point, regional, and continental-scale δ 2 H and δ 18 O values in precipitation, especially for some regions where little or no GNIP data existed.This demand for isotopic predictive capacity arose from the ecological, wildlife, food source traceability, and forensic sciences after it was revealed that the δ 2 H (and δ 18 O) values of some plant, animal, and human tissues closely mirrored isotopic patterns of precipitation (Hobson and Wassenaar, 1997;Bowen et al., 2005a).This strong water-to-biosphere isotopic linkage facilitated new areas of interdisciplinary spatial water isotope research ("isoscapes"), since δ 2 H and δ 18 O analyses could be used for determining origins of migratory species (Hobson and Wassenaar, 1997;Popa-Lisseanu et al., 2012), freshwater fish (Soto et al., 2013), food or drink products (Bowen et al., 2005b;Heaton et al., 2008), and in criminal forensics (Fraser and Meier-Augenstein, 2007;Ehleringer et al., 2008), as well as across spatial scales and regions that were hitherto problematic.For predictive isoscapes to progress further into these new areas of climatological, ecological, and hydrological research, improved and accurate spatially explicit predictive δ 2 H and δ 18   driven by and premised upon the foundational global patterns of δ 2 H and δ 18 O in precipitation, were required.
Early attempts at producing maps of δ 2 H or δ 18 O patterns in precipitation (Sheppard et al., 1969;Yurtsever and Gat, 1981) were later improved with geographical information systems (GIS) and the online availability of spatially extensive water isotope (GNIP) and climatic data sets (Rozanski et al., 1993;Dutton et al., 2005;Welker, 2012).Descriptive spatial maps of δ 2 H or δ 18 O in precipitation were first made using simple inverse distance-weighted approaches (Birks et al., 2002;Aggarwal et al., 2010).Thereafter, more complex multiple regression and interpolation prediction models were widely used (Bowen and Wilkinson, 2002;Bowen and Revenaugh, 2003).In particular, the approach of Bowen and Wilkinson -hereafter called model M1 -used fixed predictor variables (absolute latitude, squared latitude, altitude) to obtain the response variable (δ 2 H/δ 18 O in precipitation) and an interpolation parameter to optimize the fitting model.Because isotope ratios show a strong linear correlation with mean annual air temperature in nontropical regions, model M1 was able to explain 58-61 % of the isotopic variance in precipitation over a globally gridded domain.Model M1 did not downscale well particularly into tropical parts of the world, resulting in poorer fits between isotope data and model results (e.g., Africa, Asia, etc.).Some researchers therefore constrained the geographical domain (Liu et al., 2008) or added meteorological explanatory regressor variables at a global scale (van der Veer et al., 2009;Bowen, 2010).Others combined meteorological (e.g., surface air temperature, relative humidity, snowfall amount, etc.) and geospatial variables in order to obtain improved regressions/interpolations of precipitation isotope composition for specific regions such as Austria (Liebminger et al., 2006), the eastern Mediterranean (Lykoudis and Argiriou, 2007), and Canada (Delavau et al., 2011).In a different approach, a spatially dense groundwater isotope network was used as a "proxy" for mean annual precipitation and used generalized linear models (GLMs) and multiple predictor variables (elevation, precipitation amount, latitude, basin, and their interactions) to explain 81 % of the isotopic variance and patterns in precipitation in Mexico (Wassenaar et al., 2009).All of the above regionalized multivariate regression and interpolation approaches resulted in markedly improved predictive outcomes across some regions as manifested by (1) higher coefficients of determinations (R 2 ) between the measured and modeled data and (2) lower predictive uncertainties or residuals when compared to model M1 (a summary of models is given in Table 1).However, because regional models were arbitrarily fitted to small geographical domains, they were not applicable at the global scale.
The goal of this paper was to describe a new regionalized cluster-based water isotope prediction (RCWIP) approach for predictive annual, seasonal, and point location estimates of δ 2 H and δ 18 O in precipitation around the globe.In order to leverage the improved predictive accuracy accrued by regionalized multivariate approaches -and to help bridge the gap between regionalized or a one-size-fits-all global model M1 -RCWIP used fuzzy clustering of GNIP stations into predefined regionalized climatic zones.Statistical climatic zone clustering was used rather than arbitrary delineation of geospatial domains (e.g., country or geographic region).The RCWIP approach was highly flexible, because rather than applying a single fixed model to each regional domain, a suite of all possible regression models was tested, with the best performing model selected.In order to produce a spatially continuous global δ 2 H or δ 18 O isotopic map, fuzzification was used to weight and amalgamate the climatically regionalized prediction maps into a single thematic global isotopic map.A performance assessment of RCWIP was made by comparing its outputs to the established model M1, and by using the identical isotopic data sets for all models.Finally, the comparative predictions of accuracy, precision, and  uncertainty of precipitation isoscapes were used to illustrate improved RCWIP outcomes.

Materials and methods
RCWIP employed the well-established linear prediction model based on multiple regressions of station-based precipitation isotope and climatic data and the interpolation of residuals (Bowen and Revenaugh, 2003).RCWIP differed from this key earlier work in that it used the best performing model from a suite of regionalized domain multivariate regression equations to determine the isotopic composition of precipitation at a known location i as a function of the selected predictor variables available.A comparison between model M1 and the RCWIP approach is described in Fig. 2. To accommodate regionalization, RCWIP introduced additional steps of formation of subsets of station-based data as well as the gridded data sets according to their climatic clustering properties, as well as the determination of model parameter settings and the coefficient estimations for each climatic cluster, as described in turn below.

δ 2 H and δ 18 O isotope and spatial data
The precipitation stable isotope data set used was comprised of monthly composites of δ 2 H and δ 18 O in precipitation from worldwide stations archived in the GNIP database (Aggarwal et al., 2010), supplemented by compatible monthly isotope data from published papers (Wang and Peng, 2001;Kralik et al., 2003;Kurita and Ichiyanagi, 2008).The temporal span of isotopic data records was constrained to 1960-2009, although fewer than 100 stations were fully contemporaneous over this time frame.The longest serving GNIP stations spanned 20-52 yr of isotopic records.Additionally, stations that did not have over two years of monthly isotopic records were omitted to avoid seasonal biases aris-  ing from incomplete records.The calculation of amountweighted annual δ 2 H and δ 18 O composites was obtained from monthly data sets from the GNIP (IAEA/WMO, 2013).The prescreened water isotope data used comprised δ 2 H (> 49 000) and δ 18 O (> 57 000) records from GNIP stations collected between 1960 and 2009 (Fig. 1).All spatial variables (latitude, longitude, elevation) for each station were The stable isotope data in the GNIP database were composed of station-based point measurements with substantial gaps in spatio-temporal coverage (Rozanski et al., 1993;Aggarwal et al., 2010).Concerns surrounding the pooling and use of noncontiguous data sets for spatial isotope mapping have previously been discussed (Bowen, 2010).Furthermore, nearly half of the longest-serving GNIP stations showed time trends when adjusted for seasonality, although there was no obvious process to explain the trends.Additionally, for trending stations, the mean interannual rate of isotopic change was rather small (δ 18 O < 0.03 ‰ per year) in comparison to seasonal isotopic fluctuations, historical measurement error, and methodological changes, but the cumulative decadal change appeared significant over longer time frames for a few regions (Table S1 in the Supplement).For this paper, we accepted the spatio-temporal limitations of the GNIP data set in pooling noncontiguous decadal-scale data sets (Bowen, 2010).However, users of isoscape mapping efforts must remain alert to inherent assumption of temporal constancy in the isoscape maps, particularly when the studies may be sensitive to interannual or stochastic variability (e.g., plant or food traceability studies).

Climate data
Climatological variables (monthly precipitation amount, air temperature, vapor pressure) corresponding to the isotope sampling stations were also obtained from the GNIP database.When corresponding meteorological data were unavailable, averaged climate data (monthly, annual) was obtained from the nearest station in the Global Historical Climate Network station (GHCN) (Peterson et al., 1998).If a nearby GHCN station was not available, climatic parameters were obtained from the CRU-CL 2.1 data set (New et al., 2002).The precipitation amount data were natural-logtransformed to achieve normality and to linearize the relationship with δ 2 H and δ 18 O.Point-based data on monthly mean precipitable water was included into the array of climatic variables as it helped define the mean residence time of vapor, which had high predictive power for the isotopic composition of precipitation in tropical regions (Aggarwal et al., 2012).
Gridded data sets of spatial and climatic variables (mean monthly/annual temperature, precipitation, and relative humidity) were extracted from the CRU-CL 2.1 data set (New et al., 2002) at a resolution of 10 arcmin.Vapor pressure was calculated from temperature and relative humidity using the equations given by the World Meteorological Organization (WMO, 2008).Gridded precipitable water data were obtained by interpolation of point-based information; however, the use of remotely sensed data (e.g., NASA EOS Atmospheric Infrared Sounder) is anticipated in the future.

Climatic zone clustering
Because climatic parameters are spatially continuous, sharply drawn climatic classification boundaries on a map pose an unrealistic depiction of their spatially continuous distribution (McBratney and Moore, 1985).Climatic classes, however, may be useful in varied forms (e.g., plant hardiness zones).To avoid sharp boundary domain constraints, fuzzy clustering techniques, including fuzzy c-means (FCM; Bezdek, 1981;Cannon et al., 1986), were used to (1) statistically reduce the global climatic data set to a number of manageable and unique spatial domain clusters and (2) to seamlessly overlap the determined climatic class boundaries.Fuzzy clusters used fractional membership values (e.g., a GNIP station may be 60 % in cluster A and 40 % in cluster B) rather than a strict binary criterion (whether it belongs or not) to classify stations or points on the landscape.To achieve this, the FCM routine randomly seeded cluster centroids in the data space and iteratively adjusted their positions to minimize the total of all distances between the input data points and the centroids until a convergence criterion (such as minimal improvement from one iteration to the next, or preset number of iterations) was reached.Cannon et al. (1986) defined the fuzzy c-means function by with U being the fuzzy membership matrix for a set of centroids ϑ = (ϑ 1 , ϑ 2 , . .., ϑ c ).The distance d ik between the kth data point x k and the ith centroid ϑ i was calculated as (2) m ∈ [1, ∞] was a fuzzification (i.e., smoothing) factor that indicated the sharpness of transition between two fuzzy clusters, with m = 1 indicating a crisp boundary line.Cannon et al. (1986) reported a value of 1.1 < m < 5 as "useful".For RCWIP clustering, a factor of m = 1.5 was experimentally determined as a suitable compromise between fuzzy boundaries and spatial explicitness: while a smoothing factor of 1.0 would have resulted in line boundaries (no fuzzy transitions), we experimented extensively with the smoothing factor and found out that m > 1.5 would have led to inordinately large clusters extending, for example, from one continent into another ignoring the ocean in between.Since our approach was focused on land areas, we empirically found that m = 1.5 was a suitable (albeit arbitrary) compromise to balance between terrestrially meaningful clusters and obtaining smooth transitions between them.
In order to build appropriate climatic clusters, data from GHCN records (n = 5921) were used in 26-dimensional data space (26 normalized climatic and spatial variables, 12 monthly mean temperatures, precipitation, latitude, and longitude) and were clustered using weighted FCM.The number of climatic clusters was restricted to 1.5 times the number of climatic input variables, although it was generally recommended that the number of clusters be twice the number of input variables (McBratney and Moore, 1985).Our reasoning was practical and by trial and error, because we found using more climatic clusters led to difficulties in subdividing the unevenly distributed GNIP data set (Fig. 1).Each GHCN data point was then assigned a weight based on its proximity to its nearby stations, whereby a GHCN station's average distance to its five closest neighbors was normalized to this parameter's average over the whole GHCN data set.Globally, the GHCN data were heavily overrepresented in some regions, and unconstrained FCM clustering would have resulted in an excessive number of climatic clusters in weather-data-rich countries like the USA, Canada, Australia, and China (Table 2).This exercise resulted in the identification of 36 climatic clusters.
The 36 FCM-derived climatic clusters were then subjected to a first-order evaluation of their appropriateness by comparing their features to the well-known Köppen-Geiger climate classification scheme (Kottek et al., 2006).Our evaluation also consisted of determining whether each FCM-defined climate zone cluster was consistent with the geographical location of the cluster centroid.Following this, some offline manual adjustments were needed.These changes mainly involved manually moving the geographic location of FCM derived climatic centroid to the closest station of the WMO's Global Surface Network (GSN).In a few cases (e.g., Antarctica), cluster centroids were moved to entirely new positions to better cover those areas that were underrepresented by the automated FCM clustering.A spreadsheet of the FCM outputs, as well as the geographical positions and climate data averages of the 36 centroids along with descriptions of each climatic zone was tabulated in the supplementary material (Table S2).
The final positioning of the 36 global climatic zone centroids served as the fundamental basis for defining the RCWIP subsets, and was applied to the GNIP station data set as well as the gridded climatic data sets, resulting in each GNIP station or each CRU grid cell being assigned a fuzzy cluster membership u ik : (3) Membership fractions of < 0.02 were considered insignificant and omitted in order to give more weight to the main membership fractions of a data point or grid cell.The remaining fractions were renormalized to their sums.Fuzzy membership fractions were also used to layer-stack the membership of the GNIP stations to each climatic zone.The criterion used for including a GNIP station in a given climatic zone regression subset was a minimum membership fraction of 0.1.A map of the 36 climatic cluster zones used is shown in Fig. 3. Figure 4 also provides an illustrative example of the spatial depiction of clustering membership fractions and the amalgamation for a selected geographical region.This clustering procedure also provided a useful exercise for categorization of the GNIP stations, since it afforded a statistically robust and unbiased (e.g., by arbitrary boundaries) means to group subsets of "climatically similar" GNIP stations.One benefit of this analysis is that it allowed a visualization of climatically similar zones around the Earth that were (1) overrepresented by current and historical GNIP data collections (e.g., Europe), (2) lacking in GNIP data to be targeted for future stations, and (3) historically underrepresented (e.g., compare the GNIP station distribution in Fig. 1 with zones in Fig. 3).

Regression models
RCWIP tested all of the available regression-interpolation models, as summarized in Tables 1 and S3.Model M1, as noted, was the widely used global domain model that uses  (Bowen and Wilkinson, 2002;Bowen and Revenaugh, 2003).Model M1 therefore served as the basis upon which to assess RCWIP performance.Model M2 comprised the same global domain model of M1, but used a covariate best fit instead of a fixed regressor combination.To further accommodate regionalization through the use of the climatic clustering, model M3 used the same approach as M1, but using fixed regressors for each of the 36 individual climatic clusters.The most flexible model of all -model M4 -used all 36 regionalized climatic clusters and used the best fit of all regressor variable combinations for each zone.A process flowchart for RCWIP regression modeling is shown in Fig. 5. Climatically regionalized best-fit flexible regression models were derived from all possible combinations of the following variables: latitude, longitude, elevation, air temperature, precipitation amount, vapor pressure, and precipitable water.Any of the 120 combinations of these explanatory regressors were tested and ranked according to their coefficient of determination (R 2 ).For models M2 and M4, the best-fit combinations were selected if they met an acceptance criterion of R 2 of ≥ 0.5, and the number of remaining degrees of freedom was ≥ 7.5 times the number of independent variables (for M2, this was always the case).For model M4, this procedure was performed on all 36 climatic clusters; however, if no suitable regression equation satisfying the acceptance criteria was obtained for a given cluster or month, the global domain best-fit regression equation (M2) for the corresponding month was applied by default to the cluster/month.Since M2 was effectively a "backup solution" for preferred M4 model, it was not discussed further.
The "fixed-regressor" models M1 and M3 were applied to the same climatic clustering, but their equations were derived on the general form given by Eq. (4): For model M3, only those models for a given cluster month that fulfilled the same criteria as M4 were accepted, otherwise they were substituted by model M1.We found model M3 had a very high frequency of overlap with M1.Henceforth, all results and discussions hereafter concerning "regionalized" and "global" models involve only comparing the outcomes of model M4 to M1. Figure 5 depicts the decision tree for each of the M1-M4 regression models used in the RCWIP approach.
Once optimal regression models were obtained for each climatic cluster domain, these were applied to the station data set and to the gridded data, resulting in preliminary regression estimates (termed p i for the station and p x for the gridded data).For the station data set, point-based residuals were calculated using Eq. ( 5) for M1 and Eq.(6) for M4: For the latter two, the best model outcome is selected on the basis of the R 2 value.
In addition, p x of M4 was calculated as a fuzzy membership-weighted layer stack of the different regression models (Eq.7): To compare outcomes of model regression performance, models M1 (global domain, fixed regressors) and M4 (climatic regionalized, flexible regressors) were chosen for comparative illustration by illustrating the mean annual precipitation δ 18 O (δ 18 O ANN ) model outputs.

Kriging
In the final step of the RCWIP workflow, the point-based regression residuals were interpolated using Kriging and added to the regression surface.The best station-based regressions obtained were interpolated onto the gridded surface, with the Kriging results added to the regression surface.Kriging provided a robust geostatistical interpolation method exploiting the spatial autocorrelation of natural features (Matheron, 1963;Delhomme, 1978).The Kriging input function was based on the classic variogram model, which describes the degree of spatial dependence, i.e., the variance between the values within a pair of data points (x and y, cf.Eq. 8): The resulting semivariance of any particular point pair plotted against the distance between was used to create an empirical variogram, to which an experimental variogram model was fitted.The fitted function was defined by function type, nugget, sill, and range, and was used to predict values at any unknown location within a certain error range.Given the large heterogeneity of the GNIP data set (e.g., different scales for γ using different regression models, different months, stable isotopes), the recommended variogram autofitting function was augmented by a role-based definition of the modeled variogram that "visually best fit" the point clouds.Briefly, our variograms were computed with a reasonable (for hydrologic purposes) cut-off lag distance of 150 • (e.g., data points located farther away than 150 spatial degrees from each other exerted no influence on the variogram).A robust variogram estimation was used (Cressie, 1993), whereby the width of the variogram bins was set to one-hundredth of the cut-off distance.An exponential function was used to fit the point cloud.The nugget was set to the minimum of γ , and sill was defined as the minimum γ plus 2.25 times the standard deviation of γ .The range was detected automatically when the curve function was fitted to the variogram data.This parameter setting was chosen to ensure that any of the models tested (12 monthly for oxygen and hydrogen, plus one annual each, for both M1 and M4, in total 52) were subjected to the principles of identical treatment (PIT).When using the automatic parameter detection functions of the "gstat" library (Pebesma, 2004), some experimental variograms were found to be singular, with the model function not matching the curve given by the empirical var-iogram.Hence, we used a parsimonious PIT approach that would not result in singular variograms.We are well aware that this conservative PIT approach might not result in a perfect nugget/sill/range combination for each of the 52 models tested; however, comparability was maintained in all steps of the prediction/comparison workflow, which we deemed to be of greater importance for the RCWIP model comparison approach.Kriging errors were calculated as the square root of the variance of the Kriging estimator as output by the "gstat" library.Kriging model evaluations were conducted by graphical and numerical comparisons of the variograms for the two core models (model M1 vs. M4).This is a rather simple approach to quantify the comparative outcomes of different geospatial models in order to determine which gives a better result; regression uncertainty may be incorporated as additional uncertainty estimation into the future versions of the RCWIP model.

Confidence intervals
In order obtain a measure of the accuracy of the RCWIP regression and interpolation model prediction outcome, 95 % confidence intervals (CIs) were determined using the leaveone-out (n -1) "jackknife" resampling procedure (Wu, 1986;Shao and Tu, 1995).The jackknifing method determined the robustness of model results by running n i iterations, each time leaving one station data value out (i).In order to reduce computational load, these calculations were restricted to a 2 • × 2 • grid, which was sufficient for the determination of areas potentially deficient in predictive qualities.
Mathematically, the jackknife method was expressed as θJ = f (δ), where θJ stands for a jackknife estimator function of δ 18 O (or δ 2 H), and in which f (δ) encompasses both multiple regression and Kriging.A pseudo-value leaving out the station i was denoted as θ(i) = f (δ (−i) ).We calculated the mean θJ,x and standard error σ θJ,x of the jackknife for each grid cell x as the mean and standard deviation of the pseudo-values (Shao and Tu, 1995): The higher the standard error of the jackknife, the more vulnerable a certain grid cell was to leaving that station data out.This standard error used for the computation of 95 % CIs was defined as c = (±)1.96σθJ,x .

Climatic clustering of GNIP stations
RCWIP's data parameterization (i.e., minimum influence threshold, minimum ratio of remaining degrees of freedom over independent variables) required that any given climatic cluster domain needed mean values for a minimum of 15 GNIP stations for the computation of a cluster-specific regression model in model M4. Figure 6 shows the result of a histogram of the available data points per climatic cluster; the horizontal threshold line (n = 15) could be compared to the minimum number of data points required for the computation of the regionalized model M4.This plot revealed that of the 36 climatic regions defined by fuzzy clustering, 15 zones did not meet the required data threshold for applying model M4.
For these cases, model M2 was substituted by default.This graph further illustrated the relative paucity of isotopic data for specific spatial and climatic regions.For example, areas that were particularly data deficient in precipitation isotope data  USA), 33 (Pacific islands), and 34 (Australia).These and the 11 other identified areas could benefit from targeted new GNIP sampling efforts (Fig. 3).A map indicating the coverage of M4 for annual δ 18 O is available from the Supplement (S5).

Residuals analysis
Comparing the residual statistics of model M4 versus M1 in Fig. 7a, several indicators of improved performance of model M4 were identified.The standard deviation of model M4 was lower than M1 (1.58 vs. 2.31), as was the interquartile range (1.57 vs. 2.51).The excess kurtosis of M4 exceeded that of M1 (3.79 vs. 3.19).The difference between model M1 and M4 residuals ( e = e M1 − e M4 ; n = 623) revealed that > 67 % of the data point residuals were lower by using M4.A frequency distribution of e is shown graphically in Fig. M4 had lower residuals on average combined with a lower residual spread compared to M1.

Improved monthly predictions
The improved performance of model M4 in establishing regionalized precipitation isotope prediction models was particularly evident when the time domain of the GNIP data was reduced to monthly means.The coefficient of determination (R 2 ) of the regression model for any particular month was limited by the fact there were fewer underlying data.For example, model M1 generally produced fairly high R 2 values overall, but it exhibited poorer performance in some geographical areas like the tropics, and particularly when considering monthly predictions.The root-mean-squared errors (RMSE) of model M1 and M4 were examined for each calendar month and for each climatic cluster.The evaluations of all clusters revealed that the RMSE of model M4 was lower than M1 in all cases for both δ 18 O and δ 2 H, in all 36 clusters for all 12 months, and in all clusters for the mean annual model.Comparisons for several selected climatic clusters are illustrated in Fig. 8.However, given the fragmented GNIP data coverage, it should be noted that suitable model M4 regression models could only be obtained for 56 and 58 % of the δ 18 O and δ 2 H combinations due to underlying data availability constraints.
An ideal predictive isotope model would be seasonally flat and have a low RMSE.In all cases, model M4 obtained lower monthly RMSE values than model M1.In some cases, for example cluster 1 (Berlin), model M1 performed much more poorly in winter or, in the case of cluster 10 (Morocco), had the higher RMSE during the summer dry season.In short, the consistently overall lower RMSE of model M4 over M1 indicated that substantially improved seasonal and monthly predictive outcomes were obtained using the RCWIP approach.

Interpolation
Following the analysis of the multivariate regression results, an evaluation of the performance of the interpolation methods was undertaken.To achieve this, we examined the empirical and experimental semivariograms, along with the accompanying interpolation uncertainty.For brevity, only models M4 versus M1 for amount-weighted annual δ 18 O (δ 18 O ANN ) were illustrated.In and clearly indicated superior performance over model M1 (±1.38 and ±11.5 ‰, respectively).Outcomes for monthly prediction models and a table of the Kriging parameters are listed in Table S4.
Figure 10 shows a global map of the spatial distribution of the isotope prediction as the difference between predicted model outcomes of model M4 versus M1 ( σ = σ M4 −σ M1 ) for the δ 18 O ANN and δ 2 H ANN grids.The predicted isotope differences varied spatially from little to no difference between model M4 and M1 (e.g., < 0.2 ‰ for δ 18 O) to rather large differences (> 0.7 ‰ for δ 18 O).Model M4 outperformed M1 particularly in data-scarce areas (e.g., in Australia, northern and central Asia, and Greenland), but also in parts of North America and Africa.On the other hand, any advantages to using RCWIP in isotope data-rich areas like Europe appeared to be limited (usually less than 0.3 ‰ difference, but see RMSE above).

Confidence intervals
The 95 % confidence intervals (CI) for δ 18 O ANN and δ 2 H ANN models were derived as 1.96 times the standard error of the jackknifing of model M4.For δ 2 H, the spatial distribution of CIs is illustrated in Fig. 11.Note that δ 2 H is shown as an example and that δ 18 O would appear accordingly due to its strong correlation with δ 2 H (cf. Fig. 1).It was clear that 95 % CIs were below ±1 ‰ for δ 2 H (e.g., lower than measurement error) for nearly all parts of the world.However, there were regions of higher CIs in southern and western Asia, the southern parts of the Arabian Peninsula, and particularly in the vicinity of the Himalayas (e.g., mainly in climatic clusters 4, 5, 7, and 8).Higher CIs were also observed for climatic clusters 28 and 29 (southern Andes) and 32 (some Pacific islands).These higher CI anomalies may be due to the fact these climatic clusters were based on isotope data sets whose size barely exceed model M4 minimum station threshold (i.e., had there been only one or two stations less, the cluster would have deferred to model M2), or in some cases where there were extremes in geographical characteristics not well served by the spatial grid size used (e.g., very rapid elevation change in the Himalayas, high precipitation in Asia, or large distances to neighboring data points).Furthermore, CIs can also be affected workflow, whereby data-deficient areas were covered by the globally parameterized regression model M2, whose size far exceeds any of the localized climatic subsets, and thereby precludes changes in the prediction result for any given jackknifing pseudovalue.The data scarcity limitation applies, for example, to climatic clusters 16 and 18-24 (North America), 11 and 14 (western and central Africa), 34 and 35 (Australia, New Zealand), and 3 (Russia).It was anticipated that 95 % CIs may be further improved through the targeted collections of new precipitation isotope data in these regions into the future.

Global precipitation isoscape maps
The final outcome of the RCWIP approach was to construct globally gridded precipitation prediction maps (isoscape maps) of δ 18 O and δ 2 H for the Earth's land areas (excluding Antarctica), restricted to a 10 arcmin resolution.These isoscape map products included precipitation δ 18 O ANN and δ 2 H ANN , amount-weighted growing-season δ 18 O GS and δ 2 H GS , and amount-weighted monthly grids.A complete suite of isoscape map products has been made available for public use as geo-referenced TIFF (GeoTIFF) files, available from the website of the IAEA Water Resources Program (www.iaea.org/water/).These GeoTIFF files can be used in a variety of disciplinary fields of hydrologic and ecological research.
To provide a few comparative outcome examples, the predicted global isoscape map for δ 18 O ANN was determined using the RCWIP approach, and is depicted in Fig. 12.In order compare the RCWIP approach to model M1, a map showing the isotopic difference ( δ) for δ 18 O ANN and δ 2 H ANN was constructed, and is illustrated in Fig. 13.Examination of Fig. 13 reveals several substantial differences between RCWIP and model M1, in some cases for only one of the isotopes.One of the most obvious differences between RCWIP and model M1 for both isotopes was found in the Arctic regions -particularly Greenland -and for parts of northern Asia (e.g., Siberia).In these areas, RCWIP consistently produced more negative δ 2 O ANN and δ 2 H ANN predictions for precipitation.To a lesser extent, this also held true for the western Andes and for parts of northern Africa, particularly the Ethiopian Highlands (Fig. 13).Conversely, RCWIP produced more positive δ 18 O ANN and δ 2 H ANN predictions than model M1 in Australia and the Himalayas, as well as over parts of eastern Asia, to a lesser extent over southwestern North America, and in parts of the South American lowlands and the Arabian Peninsula (Fig. 13).
Decoupled isotopic differences (e.g., more enriched or depleted in 18 O than 2 H than was expected from the GMWL relationship) in the results of RCWIP compared to model M1 were observed over the Tibetan Plateau, over the Andes, and in the northern Sahara (Fig. 13).These "decoupled" isotopic differences generally occurred for climatic cluster domains for which there were sufficient isotope data to build climatecluster-specific M4 regression models for only the one but not for the other isotope (e.g., Fig. 3), especially climatic clusters 7, 8, and 29.This observation was corroborated by correspondingly higher CIs in these same areas (Fig. 11).We therefore caution isoscape map users, particularly those working in the geographical areas listed here, to carefully verify these model results with data in these areas.
Finally, there is tremendous interest in precipitation hydrogen and oxygen isoscape maps for use in forensics and plant and food authenticity, as well as in wildlife and ecological studies.Previous research showed that the δ 2 H of plant and wildlife tissue generally showed the strongest correlation with amount-weighted growing-season precipitation δ 2 H and δ 18 O as a result of seasonally relevant water uptake (growing season here defined as growth months with average temperatures > 0 • C).Thus, growing-season isoscapes were of great interest to the ecological, food authenticity, and forensics fields (Cormie et al., 1994;Hobson and Wassenaar, 1997).In Fig. 14, a global isoscape for growing-season δ 2 H is shown, and has also been provided as a GeoTIFF on the IAEA website.

Conclusions and outlook
A new regionalized cluster-based water isotope prediction (RCWIP) approach based on over 50 yr of GNIP data was demonstrated for the purposes of predicting spatio-temporal patterns of the stable isotope composition (δ 2 H, δ 18 O) of precipitation around the world.Through the use of statistically defined climatic spatial domains, a series of flexible climatic and spatial explanatory regressor variables were used for testing all available models (global-regional to fixed-flexible regressors) in order to obtain the best predictive model for each climatic cluster.The best individual cluster-optimized prediction models were seamlessly amalgamated into single global map products by using fuzzy clustering.Comparative tests revealed that the RCWIP approach outperformed the previously well-established model M1 > 67 % of the time, and furthermore clearly improved both our predictive accuracy and precision.The main outcome of the RCWIP approach was to produce improved generalized and disciplinespecific relevant and useful precipitation isoscape map products to be available for download and public use from the IAEA website.
Some precautionary notes on the use of δ 2 H and δ 18 O isoscape predictive maps must be emphasized.As noted, precipitation isoscape mapping products were based upon discontinuous long-term data sets within the GNIP database, and therefore these mapping efforts had an inherent assumption of temporal constancy in their predictive outcomes, which may not be true for some regions of the world.For many applications, small overall time changes in precipitation isotopes were not likely to be strongly manifested in some receiving environments of interest.Groundwater, for example, tended to reflect precipitation events averaged over years to decades, and so would be slowly responsive to climatically driven or stochastic weather events and isotopic changes.For other disciplines where there was a strong interest in precipitation isoscapes (e.g., plant or animal ecology, food authenticity), there may be highly relevant interan-nual, stochastic, or seasonal differences in the precipitation regimes (wet/dry years, ENSO) that could affect the amount of "relevant" precipitation water entering the soil, food webs, and biological and plant tissues (and possibly also affected the isotopic composition).For those disciplines where the timescales of "which water matters" was critical, the use of isoscape map products may be useful as a first-order starting point for predictive modeling, but should be used with great caution and need to be tested and validated (e.g., annual calibrations with GNIP data).
Regarding spatial and temporal data coverage, RCWIP provided a strong and a flexible platform for predicting precipitation isoscapes, although modeling efforts were only as good as the supporting foundational data pillars.Spatial gaps in precipitation isotope collected data at global scales, spatially and temporally, over many decades were inevitable.Here, the climatic fuzzy clustering of GNIP stations allowed us to identify particularly relevant and data-deficient areas of the world.RCWIP thereby provided a basis upon which to better inform future volunteer efforts in contributing to the GNIP database.In particular, new multiyear efforts of strategically located GNIP collections would be useful in certain areas of Africa, central Asia, and in South America.
Although the RCWIP mapping effort presented here was temporally restricted from the 1960s to 2009, current and upto-date isoscape maps will be published online on the IAEA Water Resources website, with regular updating of RCWIP to improve both the global spatial coverage of GNIP and predictive modeling performance outcomes into the future.This will facilitate improved isoscape modeling for a host of disciplines and new research.Finally, we encourage users of isoscape products to become engaged in volunteer efforts to improve the GNIP coverage, both spatially and temporally, particularly in those areas that were identified to be data deficient.
O models, Published by Copernicus Publications on behalf of the European Geosciences Union.Paper | Discussion Paper | Discussion Paper | Discussion Paper |

a
Spatial variables shall be defined as any combination of the following: latitude, longitude, and elevation (including derived parameters like squared or absolute latitude).b climatic variables are defined as any combination of the following: precipitation amount, air temperature, vapor pressure and other climate-related variables (wind speed, snow amount, relative humidity, precipitable water/moisture residence time, etc.).

Fig. 2 .
Fig. 2. Workflow of the regression-interpolation models for predicting the isotopic composition of p cipitation around the globe.The upper panel (a) shows the workflow of the fixed regressor, global dom Model M1 (Bowen and Wilkinson, 2003).The lower panel (b) shows the workflow of RCWIP with addition of climate data, regionalization and the use of flexible regressors (Models M2-M4).

Fig. 2 .
Fig. 2. Workflow of the regression-interpolation models for predicting the isotopic composition of precipitation around the globe.The upper panel (a) shows the workflow of the fixed regressor, global domain model M1 (Bowen and Wilkinson, 2002).The lower panel (b) shows the workflow of RCWIP with the addition of climate data, regionalization, and the use of flexible regressors (models M2-M4).

Fig. 3 .
Fig. 3. Map of 36 climatic zone domains used in RCWIP, derived by fuzzy clustering but here shown with defined boundaries.The "+" symbol denotes the spatial location of each climatic zone centroid.The descriptions of each of these climatic zones are denoted by a number (Zone #) and a Köppen-Geiger climatic zone description (e.g.Dfc).These are fully tabulated inTable 2 and Sheet S2 in the Supplement.

Fig. 4 .
Fig. 4. Example of fuzzy climatic zone layering and amalgamation for a geographical domain around the Mediterranean-Black Sea region.The five climatic clusters in this spatial domain (climatic clusters 1, 2, 4, 10, and 12) were overlain, with degree of color saturation of each indicative of the fractional membership for each grid cell.The GNIP station in Athens (indicated by a small star sign), Greece, for example, holds membership fractions of 0.14, 0.34, and 0.52 in clusters 1, 10, and 12, respectively.

Fig. 5 . 34 Fig. 5 .
Fig. 5. Flowchart of available regression models built into RCWIP.The left portion of the diag the fixed explanatory variable regressor Models M1 and M3, using a global or climatic clus The right portion of the diagram depicts the flexible regressor Models M2 and M4 using climatic cluster domains.For the latter two, the best model outcome is selected on the bas value.

S.Fig. 6 .Fig. 6 .
Fig. 6.GNIP station data points per climatic cluster.The blue line indicates the minimum number of data points (15) to successfully derive a M4 best-fit regression model for a given cluster.
(up to 2009)  to allow application of model M4 were climatic clusters 19 (western coast of Canada), 22 (western

Fig. 8 .
Fig. 8. Regression model (δ 18 O) performance comparison based on RMSE (‰, left axis) of models M1 (red) versus M4 (blue) for amountweighted monthly predictions for six climatic cluster centroids.Lower RMSE means better performance.Mean monthly precipitation amount (mm month −1 , right axis) is shown by blue bars to indicate the timing of rainy and dry seasons.An ideal model would produce seasonally consistently low RMSE.That there were flatter monthly patterns and a lower RMSE of M4 compared to M1 shows improved seasonal predictive outcomes when using regionalized best-fit regressor combinations.

Fig. 11 .Fig. 12 .
Fig. 11.The 95 % CIs of Model M4 for mean δ 2 H ANN in precipitation.The legend scale is ±δ 2 H expressed in ‰.The highest CIs (poorest performance) was clearly observed in the Himalayan region in Asia.See text for discussion.

Fig. 13 .
Fig. 13.Comparison of difference ( δ = δ M4 − δ M1 ) in predicted amount-weighted δ 18 O ANN (upper panel, scale) and δ 2 H ANN (lower panel, scale) comparing model M1 and M4.The isotopic scale bar and color scheme indicate the extent -positive or negative -and locations where M4 produced more positive or negative predicted values than M1.

S.
Terzer et al.: Global isoscapes for δ 18 O and δ 2 H in precipitation aper | Discussion Paper | Discussion Paper | Discussion Paper |

Fig. 14 .Fig. 14 .
Fig.14.Predicted amount-weighted mean growing season δ 2 H GS in precipitation obtained using RCWIP.Growing season was defined as the mean of all months where the average air temperature was > 0 • C. Legend is δ 2 H in ‰, relative to VSMOW-SLAP scales.

Table 1 .
Examples of global and regional water isotope models and type of predictive spatial and climatic variables used for the isotopic composition of precipitation and ground water.

Table 2 .
The ID number, geographic centroid location, climatic zone, and WMO station of the climatic clusters used in RCWIP.