How reliable are modeled precipitation isoscapes over a high-relief mountainous region?

Precipitation isotope maps over mountains are critical for water resource assessment, yet isoscape modeling has been minimally investigated for complex terrain with high relief. Here we show that multiannual (2011–2015) mean precipitation isoscapes across the Japanese Alps region can be represented by simple multiple-regression models with strong goodness of fit (R2 = 0.928 for δ2H and 0.944 for δ18O). Reliability of the models was further confirmed by agreement with previously reported independent data throughout a wider range of elevation (7–3,730 m asl). Modeled precipitation isoscapes were consistent with observations of soil water (20 sites) and river water (50 sites) when considering evaporative enrichment and unclosed water balance. Uncertainties of modeled δ values for precipitation were greater in lowlands near the coast than over the mountains. Unexpectedly, spatially uniform isotopic lapse rates (–11.66‰ km–1 for δ2H and –1.724‰ km–1 for δ18O), which are divorced from continental effects, gave substantially good approximation on a regional (e.g., a few hundred kilometers) scale. This allows for convenient and reliable modeling of precipitation isoscapes over mountainous regions for hydrological/ecological/interdisciplinary applications.


INTRODUCTION
Stable isotopes represent a powerful tool for tracing flow/ migration pathways and formation/modification processes of various substances over a wide range of spatial scales (e.g., global, continental, catchment).Isotope maps often called "isoscapes" are essential for many applications in atmospheric science, hydrology, geology, plant and animal ecology, microbiology, anthropology, and forensic sciences (West et al., 2010;Bowen, 2010a).Among others, stable isotopes of oxygen and hydrogen in meteoric precipitation are the most fundamental and widely applied (Dansgaard, 1964;Bowen, 2010a).
A major remaining issue is reliability of isoscape models for high-relief mountainous regions.It is known that precipitation isotope ratios are affected by topography (especially elevation), and isotopic lapse rate (i.e., depression of isotopic δ values per unit increase in elevation) varies with location (Clark and Fritz, 1997).Therefore, the spatial distribution of precipitation isotope ratios over mountainous regions is expected to be very complex, reflecting topography with strong relief.Indeed, West et al. (2014) reported large discrepancies between observed and modeled isotope ratios in high elevation regions of South Africa.Improvements in precipitation isoscape modeling over mountains as natural "water towers" (Viviroli et al., 2007) are critical for water resource assessment.Nevertheless, there have been few studies of precipitation isoscapes focused on mountainous regions, except for the Alps (Liebminger et al., 2006a, b;Kern et al., 2014).
Another issue is the evaluation of uncertainty in values predicted by isoscape models.In statistical models, although this uncertainty is usually evaluated by goodness-of-fit (to calibration data), few studies have validated the models using data independent of model development (Lykoudis and Argiriou, 2007).Geostatistical models do not allow for any estimates of prediction errors because they are optimized to completely fit modeled isotope ratios to observed ones (Bowen, 2010b).Thus, the validity of interpolated values between observation sites has not been sufficiently examined.This problem may be serious, especially in mountainous areas, because strong topographic relief may introduce extreme local change in precipitation isotope ratios.
The objectives of the present study are to develop a precipitation isoscape model for the highest-relief mountainous region in Japan, and confirm its reliability/uncertainty from the following three aspects: (1) goodness of fit; (2) agreement of interpolated/extrapolated values with independent observation data; and (3) consistency of the precipitation isoscape with soil/river water isoscapes.Here we focus on multiannual averages of precipitation-amount-weighted mean values.This is because these are essential in hydrologic applications assuming steady-state conditions, such as estimation of groundwater recharge elevations (Yasuhara et al., 1993) or areas (Tsujimura et al., 2007), mixing analysis of various recharge sources (Liu and Yamanaka, 2012), evaluation of lake water balance (Gibson, 2001), and partitioning of local and nonlocal water resource use (Good et al., 2014).

MATERIALS AND METHODS
For a 5-year period (or shorter at some sites but always longer than 2 years; see Supplement Table SI) from July 2011 through June 2015, precipitation isotopes were monitored at 16 sites (Figure 1) across a mountainous region of central Japan, including the Japanese Alps and many volca-noes (e.g., Mt.Fuji, Mt.Ontake, Mt.Yatsugatake).Elevations of the sampling sites were 304-1,594 m asl (above sea level).Precipitation samples were collected quasi-monthly (i.e., not identical with calendar months, and often variable by site) using precipitation collectors specially developed to avoid evaporation of stored water (Yamanaka et al., 2005) using a polyethylene funnel (8.5 cm in diameter).From December through March, the funnel was replaced by a metal one (7.0cm in diameter), whose temperature was maintained at 3°C by electronic heating to melt snow.Although measurements of precipitation amount at mountain sites are affected by strong wind, no such an effect has been reported for isotopes and thus was neglected.
In addition to precipitation sampling, we sampled soil water at 20 sites (Supplement Table SII) within the Chikuma River (upper Shinano River) and Fuji River basins, and river water at 50 sites (Supplement Table SIII) for those rivers and their tributaries.River water sampling began in August 2010 at 3-month intervals for 2 years (1 year for soil water).Elevations of the sampling sites were 6-2,327 m asl for soil water and 19-2,121 m asl for river water.Soil water samples were extracted by centrifugation or cryogenic vacuum distillation methods from soil cores of 1 cm in diameter and 90 cm in depth from the ground surface.Soil cores were collected by an auger at forested sites.
Hydrogen and oxygen stable isotope ratios ( 2 H/ 1 H and 18 O/ 16 O) in all samples were measured by tunable-diode laser spectroscopy with a liquid water isotope analyzer (L1102-i, Picarro Inc., Santa Clara, CA, USA).Measurement results are expressed using common δ notation (δ 2 H and δ 18 O).Measurement errors for the analyzer were 1.0‰ for δ 2 H and 0.1‰ for δ 18 O (Yamanaka and Onda, 2011).A preliminary report on the first 1.5 years of monitoring focusing on deuterium excess d (≡ δ 2 H-8δ 18 O) of precipitation is detailed in Wakiyama et al. (2013).
For the precipitation isotope ratio (δ PPT ), we calculated precipitation-amount-weighted 12-month moving averages for each site and then obtained a long-term (multiannual) mean value by averaging them to minimize errors caused by seasonality.For river and soil water isotope ratios, arithmetic means were used as long-term means.Details of data processing and analysis are given in Supplement Text S1.
For modeling precipitation isoscapes, the following were used as potential explanatory variables in multiple regression analysis (MRA): latitude, latitude squared, longitude, longitude squared, elevation, slope, slope aspect, slope curvature, annual precipitation, annual mean temperature, sunshine duration, annual solar radiation, and annual maximum snow depth.The reason latitude squared and longitude squared were included is the latitudinal or longitudinal dependence of δ may not be always linear.Additionally, we supposed that sunshine duration or annual solar radiation may reflect precipitation frequency and its seasonality.Stepwise MRA was adopted to select essential variables with a criterion of p < 0.05 for either δ 2 H or δ 18 O.In advance of MRA annual mean temperature was excluded out of potential explanatory variables to avoid multicollinearity with elevation.
As independent data for validating modeled precipitation isoscapes, we gathered and compiled data of annual weightedmean precipitation isotope ratio reported in the literature (33 sites in total; Supplement Table SIV).The selected sites cover wider ranges of elevation (7-3,730 m asl), latitude Figure 1.Map of the study region and sampling sites of precipitation (PPT), soil water (SW) and river water (RW) (34.8°-37.1°N)and longitude (136.2°-139.4°E)than those of our original sites.
Although measurement errors are 1‰ for δ 2 H and 0.1‰ for δ 18 O, these values are for a single sample.The accuracy of long-term means should be higher, and decimal fractions affect coefficients of MRA-derived equations.Thus, we report δ values with a greater number of decimal digits than the error for a single sample.
Elevation dependence of the multiannual mean δ PPT was significant with correlation coefficients r = 0.750 (p < 0.001) for δ 2 H and r = 0.806 (p < 0.001) for δ 18 O.Stepwise MRAderived hydrogen and oxygen isoscape models for precipitation were expressed using the following functions: The δ PPT values modeled by Equations ( 1) and ( 2) were in good agreement with those independently observed in other studies (Figure 2), with r = 0.925 (p < 0.001) for δ 2 H and 0.935 for δ 18 O (p < 0.001).The MAE (or RMSE) was 4.77‰ (5.72‰) for δ 2 H and 0.631‰ (0.865‰) for δ 18 O.These are two to four times greater than those of Equations ( 1) and (2).Since the observation period was limited to only one year in most previous studies, the above measures of uncertainty include both errors in the regression and those from interannual variability in observed δ PPT .Weak negative (i.e., underestimation) bias in modeled δ 2 H was found in northern areas along the coast of the Sea of Japan.Particularly large deviations in both modeled δ 2 H (> 15‰) and δ 18 O (> 1.5‰) were evident at lowland sites near the Pacific coast, on the southwestern and southeastern edges of the study region (Figure 3).Although this validation includes extrapolated values (i.e., the elevation range exceeds that for our own data used in MRA), Equations ( 1) and (2) work well.For instance, errors at the highest site (3,730 m asl) near the summit of Mt.Fuji were only 4.40‰ for δ 2 H and 0.038‰ for δ 18 O.
Similarly, modeled δ PPT averaged over catchments corresponding to a river water sampling site tended to be lower than observed isotope ratios for river water, δ RW (Supplement Figure S2), with R 2 of 0.809 for δ 2 H and 0.715 for δ 18 O.Mean bias was -4.10‰ for δ 2 H and -0.429‰ for δ 18 O, with a δ 2 H -δ 18 O slope of 9.6.

Isotopic lapse rate
Although the spatial pattern of the precipitation isoscape was very complex (Supplement Figure S3), it could be mod- eled by simple functions of elevation, latitude and longitude, and the elevation dependence (i.e., "altitude effect") was among the most significant.In other words, the complexity in precipitation isoscape is attributable to complex terrain, and rules governing the isoscapes are simple.Accuracy of the modeled δ PPT did not depend on elevation and was reasonably good even in alpine zones near summits.These results suggest that a spatially uniform isotopic lapse rate is a good approximation for a high-relief mountainous region on a horizontal scale of a few hundred kilometers.
Values of the isotopic lapse rate (-11.66‰km -1 for δ 2 H and -1.724‰ km -1 for δ 18 O) are nearly the same as, or somewhat smaller than, those reported in other studies (e.g., Poage and Chamberlain, 2001;Supplement Table SV).In general, the Rayleigh process is a fundamental cause of not only the altitude effect but also the continental and latitude effects (Clark and Fritz, 1997), such that these isotopic effects are difficult to fully distinguish (Ingraham, 1998).The latitude dependence of δ PPT including a quadratic term in Equations ( 1) and ( 2) appears to be a reflection of the continental effect, because modeled δ PPT for a fixed elevation and longitude has a minimum at latitude ~36.3°N, distant from both the coasts of the Pacific and Sea of Japan.Because air masses producing precipitation generally move eastward because of the prevailing westerlies, the longitude depen-dence of modeled δ PPT (i.e., decrease with increasing longitude) can be a part of the continental effect.The present study evaluated the altitude and continental effects separately, while most previous studies measured isotopic lapse rates that were not divorced from continental or other isotopic effects.That may be the reason why isotopic lapse rates in the present study are smaller than those in the literature.This is most likely the case in relatively warm climates, although the rates may be actually larger in colder climates (Gonfiantini et al., 2001).
According to our preliminary assessment on isoscape modeling using the first 1.5 years of monitoring data (Makino, 2013), topographic and climatic variables (slope, slope aspect, annual precipitation, and maximum snow depth) were required to effectively model precipitation isoscapes by MRA.However, the present study, using longer-term (i.e., 2-5 years) data, derived simpler but more accurate models.As for the European Alps, annual or multiannual mean precipitation isoscapes showed strong coherence with elevation, whereas seasonal isoscapes require more complex functions (including climatic variables) to be modeled (Liebminger et al., 2006a, b;Kern et al., 2014).These facts imply that the effects of local and/or transient processes on long-term precipitation isoscapes cancel each other, and the altitude effect with constant lapse rate becomes more significant.Regionally constant lapse rates are evident in observed isoscapes of stream water in mountainous region (Brooks et al., 2012;Katsuyama et al., 2015), probably for the same reason as above; catchment transit time can exceed several years (Ma and Yamanaka, 2013).

Soil and river water isoscapes
Reliability of modeled δ PPT is supported, in part, by consistency with δ SW and δ RW , although there were non-negligible biases.In general, soil water and shallow groundwater can be more enriched or depleted in heavier isotopes than annual weighted mean precipitation.Isotopic enrichment is primarily caused by evaporation of water before percolation, and depletion is mainly a reflection of covariance in the seasonality of recharge flux and isotopic composition of precipitation (Clark and Fritz, 1997).
For a forested site near Mt.Yatsugatake in the present study region, Tsujimura and Tanaka (1998) reported that soil water at depths 0-1 m had δ 2 H (δ 18 O) values 2.2-3.7‰(0.39-0.57‰) larger than that of the weighted annual mean for precipitation, and attributed this enrichment to isotopic fractionation caused by soil surface evaporation.These values of enrichment are comparable with our results.In addition, the δ 2 H-δ 18 O slope for a mean bias of 5.2 is equivalent to that produced by kinetic fractionation caused by evaporation under a relative humidity of 75% (Gonfiantini, 1986;Clark and Fritz, 1997), which is reasonable for this region with a temperate humid climate.These facts indicate that the difference between observed δ SW and modeled δ PPT is not just a systematic error but an actual effect of isotopic fractionation due to evaporation.
Factors controlling the isotope separation of river water from precipitation, Δ RW (defined as observed δ RW minus catchment mean modeled δ PPT ), are potentially different from those for soil water, because the δ 2 H-δ 18 O slope of 9.6 is too large.Figure 4a shows strong correlation between observed deuterium excess for river water (d RW ) and the areal fraction of paddy fields and open water bodies (e.g., rivers, lakes, reservoirs) within each catchment (F P+W ).This was so, except for northern catchments where maximum snow depth exceeded 0.5 m, suggesting evaporative enrichment.In Figure 4b, although the Δ RW values are large in catchments with large F P+W , they are not always smaller with smaller F P+W .This indicates that there are mechanisms or reasons for increasing Δ RW , together with evaporative enrichment.
Although one potential reason is inappropriate modeling of δ PPT , it is not likely to introduce a large bias of >10‰ on catchment average.Another cause may be inappropriate modeling of the stream network, which is the basis for catchment delineation.As far as we could confirm (see Supplement Text S1), however, there were only minor errors, which cannot explain the large bias.The other possibility is an unclosed water balance or groundwater outflow across the divide.We assumed that all water precipitated onto a catchment runs off as rivers/streams at the sampling site, although some portion may run off from catchments as groundwater.If waters recharged at higher (rather than lower) elevations tend to run off as groundwater through deep flow paths (as illustrated by Tóth, 1963), then modeled δ RW is underestimated and Δ RW increases as a result.This would be the case in some small catchments, because regional groundwater flow across the divide of such catchments tends to be significant in high relief terrain (Gleeson and Manning, 2008).Although other factors such as temporal or spatial coherence between δ PPT and recharge flux potentially affect δ RW (Bowen et al., 2011), they do not always increase or tend to decrease Δ RW ; the negative Δ RW at a site RC6 (see Supplement Table SIII) may be explained by this effect.Thus, it is likely that a larger Δ RW is mainly attributed to an unclosed water balance as well as evaporative enrichment.

CONCLUSIONS
We developed simple, multiple regression models of longterm mean precipitation isoscapes over a mountainous region.Reliability of the models was confirmed by not only strong goodness of fit to calibration data but also substantial agreement with independent validation data, with MAEs two to several times as large as measurement errors.In addition, modeled precipitation isoscapes were consistent with observed isoscapes of soil and river waters, if we consider the potential effects of evaporative enrichment and an unclosed water balance within a catchment.Uncertainties of modeled δ PPT tended to be greater in lowlands near the coast.Those limitations are attributed to incomplete representation of the continental effect rather than the altitude effect.
The most important finding is that spatially uniform isotopic lapse rates give substantially good approximation on a regional scale, despite the complex terrain with high relief.This ensures convenient and reliable modeling of precipitation isoscapes over mountainous regions without geostatistical estimation, and makes isoscape approaches more useful in many applications to such regions.

Figure 3 .
Figure 3. Spatial distribution of errors in modeled a) δ 2 H and b) δ 18 O for precipitation.Negative errors mean underestimations.Stars and circles represent data from this study and previous studies, respectively

Figure 4 .
Figure 4. Relationship of areal fraction of paddy fields and open water bodies (F P+W ) with a) observed deuterium excess for river water (d RW ) and b) isotope separation for river water (Δ 2 H RW and Δ 18 O RW ; see text for definition).A regression line in a) and data plots in b) are for catchments with catchment mean values of annual maximum snow depth (Z ms ) ≤ 0.5 m

Table SI .
Summary of observed precipitation isotope data Table SII.Summary of observed soil-water isotope data Table SIII.Summary of observed river-water isotope data Table SIV.Summary of precipitation isotope data previously reported in the literature

Table SV .
Summary of isotopic lapse rates for precipitation (not including snowpack, firm, ice and surface/subsurface waters) reported in the literature FigureS1.Comparison of modeled δ for precipitation (δ PPT )