Analysis of the Spatial Distribution of Stable Oxygen and Hydrogen Isotopes in Precipitation across the Iberian Peninsula

: The isotopic composition of precipitation provides insight into the origin of water vapor, and the conditions attained during condensation and precipitation. Thus, the spatial variation of oxygen and hydrogen stable isotope composition ( δ p ) and d-excess of precipitation was explored across the Iberian Peninsula for October 2002–September 2003 with 24 monitoring stations of the Global Network of Isotopes in Precipitation (GNIP), and for October 2004–June 2006, in which 13 GNIP stations were merged with 21 monitoring stations from a regional network in NW Iberia. Spatial autocorrelation structure of monthly and amount weighted seasonal / annual mean δ p values was modelled, and two isoscapes were derived for stable oxygen and hydrogen isotopes in precipitation with regression kriging. Only using the GNIP sampling network, no spatial autocorrelation structure of δ p could have been determined due to the scarcity of the network. However, in the case of the merged GNIP and NW dataset, for δ p a spatial sampling range of ~450 km in planar distance (corresponding to ~340 km in geodetic distance) was determined. The range of δ p , which also broadly corresponds to the range of the d-excess, probably refers to the spatially variable moisture contribution of the western, Atlantic-dominated, and eastern, Mediterranean-dominated domain of the Iberian Peninsula. The estimation error of the presented Iberian precipitation isoscapes, both for oxygen and hydrogen, is smaller than the ones that were reported for the regional subset of one of the most widely used global model, suggesting that the current regional model provides a higher predictive power.


Introduction
Since the discovery that isotopic composition of precipitation can give an insight into the origin of water vapor, the conditions attained during condensation and precipitation [1], stable water isotopes have become important natural tracers in the study of the water cycle [2,3]. Isotope ratios in precipitation mainly depend on geographic position and climate conditions [1,4]. In Europe, the continental, altitude, and latitude effects are the most important [5].
Our understanding of the behavior of water isotopes has rapidly improved thanks to recent advances in the measurement and modeling of stable water isotopes [6]. The spatial autocorrelation Figure 1. The Iberian Peninsula, with precipitation stable isotope measurement networks and Köppen climate zones [34]. The 24 stations of the Global Network of Isotopes in Precipitation (GNIP) network [23] functioning between October 2002-October 2007 (red solid circles), and two regional networks (NW: [35]; and SE: [36]) are marked with blue solid circles and black circles, respectively (A). The seasonal cycle of δ 18 Op and d-excess of precipitation of the four characteristic climate zones (B, C, D, E) was computed for mean monthly data of a common six-yr-long (2001)(2002)(2003)(2004)(2005)(2006) reference period. The color coded Köppen climate zones [34] overlaid on the World_Terrain_Base basemap (Sources: Esri, USGS, NOAA) in ArcMap 10.5 [37]. The red area on the Europe inset map shows the Iberian Peninsula.
The opportunity arose to expand the GNIP dataset with precipitation stable isotope data from two regional monitoring campaigns ( Figures 1A and 2). In the first case, 21 stations were available in NW Iberia (October 2004-June 2006 [35]. In the second case, 10 stations (May 2006) and 19 stations (April 2007) were available from SE Iberia [36]. The SE regional network is located in a spatially more confined mountainous area with abrupt differences in elevation between neighboring sites, while that in the NW is better distributed, both vertically and horizontally over a larger region. Unfortunately, merging the GNIP and SE Iberia datasets drew attention to a critical lack of spatial homogeneity ( Figures S2-S4), and this turned out to be unsuitable for use in the geostatistical analysis (for details, see Electronic Supplementary Material), unlike the combination of the GNIP and the NW network. In the latter case, the average number of stations with available records from the GNIP database (~13) increased to ~34 for the period October 2004 to June 2006 (Figures 1 and 2). The Iberian Peninsula, with precipitation stable isotope measurement networks and Köppen climate zones [34]. The 24 stations of the Global Network of Isotopes in Precipitation (GNIP) network [23] functioning between October 2002-October 2007 (red solid circles), and two regional networks (NW: [35]; and SE: [36]) are marked with blue solid circles and black circles, respectively (A). The seasonal cycle of δ 18 O p and d-excess of precipitation of the four characteristic climate zones (B, C, D, E) was computed for mean monthly data of a common six-yr-long (2001)(2002)(2003)(2004)(2005)(2006) reference period. The color coded Köppen climate zones [34] overlaid on the World_Terrain_Base basemap (Sources: Esri, USGS, NOAA) in ArcMap 10.5 [37]. The red area on the Europe inset map shows the Iberian Peninsula.
The opportunity arose to expand the GNIP dataset with precipitation stable isotope data from two regional monitoring campaigns ( Figures 1A and 2). In the first case, 21 stations were available in NW Iberia (October 2004-June 2006) [35]. In the second case, 10 stations (May 2006) and 19 stations (April 2007) were available from SE Iberia [36]. The SE regional network is located in a spatially more confined mountainous area with abrupt differences in elevation between neighboring sites, while that in the NW is better distributed, both vertically and horizontally over a larger region. Unfortunately, merging the GNIP and SE Iberia datasets drew attention to a critical lack of spatial homogeneity ( Figures S2-S4), and this turned out to be unsuitable for use in the geostatistical analysis (for details, see Electronic Supplementary Materials), unlike the combination of the GNIP and the NW network. In the latter case, the average number of stations with available records from the GNIP database (~13) increased to~34 for the period October 2004 to June 2006 (Figures 1 and 2).
The precipitation records archived alongside monthly δ p values in the GNIP database were found in some cases to have gaps and/or be erroneous, raising questions regarding the calculation of seasonal and annual precipitation-amount weighted averages from monthly δ p values. Therefore, gridded (0.5 • × 0.5 • ) monthly rainfall totals from the Global Precipitation Climatology Centre (GPCC) database [38], derived as precipitation anomalies at stations interpolated and then superimposed on the GPCC Climatology V2011 [39], were utilized instead. Nevertheless, there were some discrepancies between the monthly precipitation in the GNIP records and the corresponding GPCC values: r min = 0.82 (Gibraltar), r max = 0.99 (Almeria Aeropuerto), and r mean = 0.95 (17 GNIP stations) for 2000-2013. An obvious advantage of the gridded data set is the lack of gaps and the fact that it is less prone to single-station error. The precipitation records archived alongside monthly δp values in the GNIP database were found in some cases to have gaps and/or be erroneous, raising questions regarding the calculation of seasonal and annual precipitation-amount weighted averages from monthly δp values. Therefore, gridded (0.5° × 0.5°) monthly rainfall totals from the Global Precipitation Climatology Centre (GPCC) database [38], derived as precipitation anomalies at stations interpolated and then superimposed on the GPCC Climatology V2011 [39], were utilized instead. Nevertheless, there were some discrepancies between the monthly precipitation in the GNIP records and the corresponding GPCC values: rmin = 0.82 (Gibraltar), rmax = 0.99 (Almeria Aeropuerto), and rmean = 0.95 (17 GNIP stations) for 2000-2013. An obvious advantage of the gridded data set is the lack of gaps and the fact that it is less prone to single-station error.

Steps of the Analysis and Preprocessing
The first step in data preprocessing was to numerically check the δp values for database errors (e.g. typos, sign errors) [40]. For example, very negative d-excess values, occasionally even lower than −10‰, corresponding to a monthly precipitation total exceeding 5 mm was all carefully compared to the surrounding stations. The negative d-excess value was interpreted as evidence of evaporative enrichment of the sample and it was discarded from the evaluation if no neighboring stations reported similarly extreme d-excess value for the given month. In addition, in July 2003 at GNIP station Madrid, somewhat improbably, the same values were found as in June, without any precipitation in July being recorded at that site. Consequently, the fact that the June values reappeared in July was suspected to be a database error, and was thus omitted.  [35] in the latter period. The purple lines indicate the number of records of the merged GNIP and SE network [36] for May 2006 and Apr 2007, resulting in 24 and 33 records, respectively.

Steps of the Analysis and Preprocessing
The first step in data preprocessing was to numerically check the δ p values for database errors (e.g., typos, sign errors) [40]. For example, very negative d-excess values, occasionally even lower than −10% , corresponding to a monthly precipitation total exceeding 5 mm was all carefully compared to the surrounding stations. The negative d-excess value was interpreted as evidence of evaporative enrichment of the sample and it was discarded from the evaluation if no neighboring stations reported similarly extreme d-excess value for the given month. In addition, in July 2003 at GNIP station Madrid, somewhat improbably, the same values were found as in June, without any precipitation in July being recorded at that site. Consequently, the fact that the June values reappeared in July was suspected to be a database error, and was thus omitted.
Due to the low amount, or frequently the complete lack of precipitation in summer in southern and eastern Iberia (Figure 1), the May to August monthly data were discarded from the analysis (grey bars in Figure 3). Nevertheless, the summer (defined as April to September) seasonal average was calculated, despite the unrealistically low values of d-excess and high values of δ p , because the corresponding precipitation amounts were also low, having only a small impact on the amount weighted summer seasonal value.
Due to the low amount, or frequently the complete lack of precipitation in summer in southern and eastern Iberia (Figure 1), the May to August monthly data were discarded from the analysis (grey bars in Figure 3). Nevertheless, the summer (defined as April to September) seasonal average was calculated, despite the unrealistically low values of d-excess and high values of δp, because the corresponding precipitation amounts were also low, having only a small impact on the amount weighted summer seasonal value. An amount-weighted mean was only calculated for a period (season, year) if at least 75% of the fallen precipitation was analyzed for the given isotope. This required completeness is a bit stricter criterion than the GNIP protocol (70%) [40].

Multiple Regression Analysis
Trend-like tendencies in δp over large distances across the Iberian Peninsula have been documented [24,29,35]. These empirical relationships of δp with geographical factors (e.g. altitude effect) may mask the finer scale spatial autocorrelation patterns, which are the focus of interest in the present investigation. The effect of the geographical factors influencing the water stable isotopes' variability in precipitation first has to be minimized in order to obtain representative results with variography (Section 2.3.3), i.e. the spatial trend has to be removed from the data to obtain secondorder stationarity [41]. Thus, the effect that was linked to the geographical factors was minimized by computing best-fit multiple regression models and determining their residuals.
The investigation was conducted while using the stepwise procedure suggested by O'Brien, with latitude, longitude (Web Mercator EPSG: 3857), elevation, and geodetic distance from the coast as the independent predictor variables, and the primary precipitation stable isotope values as the dependent variables [42]. Of the four independent variables, the first three were found to be applicable to the minimalization of the effect of geographical factors influencing the raw δp records. Probabilities ranged between p < 2.9 × 10 and p < 0.04, variance inflation was negligible (Variance Inflation Factor < 1.187), and the average adjusted values of R 2 were 0.63 and 0.58 for δ 18 Op and δ 2 Hp, respectively (Table S1). No trend removal was found to be necessary in the case of the derived isotopic parameter (d-excess). A similar multiple regression approach was applied to determine the degree of The number of precipitation stable isotope records in the initial monthly and multi-monthly datasets before and after preprocessing entered into the semivariogram analysis. Winter: October-March; summer: April-September. Grey bars represent the months omitted from the analysis, see text for details. The open part of a column indicates that some values were discarded after the quality check of the raw data (see Section 2.3.2).
An amount-weighted mean was only calculated for a period (season, year) if at least 75% of the fallen precipitation was analyzed for the given isotope. This required completeness is a bit stricter criterion than the GNIP protocol (70%) [40].

Multiple Regression Analysis
Trend-like tendencies in δ p over large distances across the Iberian Peninsula have been documented [24,29,35]. These empirical relationships of δ p with geographical factors (e.g., altitude effect) may mask the finer scale spatial autocorrelation patterns, which are the focus of interest in the present investigation. The effect of the geographical factors influencing the water stable isotopes' variability in precipitation first has to be minimized in order to obtain representative results with variography (Section 2.3.3), i.e. the spatial trend has to be removed from the data to obtain second-order stationarity [41]. Thus, the effect that was linked to the geographical factors was minimized by computing best-fit multiple regression models and determining their residuals.
The investigation was conducted while using the stepwise procedure suggested by O'Brien, with latitude, longitude (Web Mercator EPSG: 3857), elevation, and geodetic distance from the coast as the independent predictor variables, and the primary precipitation stable isotope values as the dependent variables [42]. Of the four independent variables, the first three were found to be applicable to the minimalization of the effect of geographical factors influencing the raw δ p records. Probabilities ranged between p < 2.9 × 10 −8 and p < 0.04, variance inflation was negligible (Variance Inflation Factor < 1.187), and the average adjusted values of R 2 were 0.63 and 0.58 for δ 18 O p and δ 2 H p , respectively (Table S1). No trend removal was found to be necessary in the case of the derived isotopic parameter (d-excess). A similar multiple regression approach was applied to determine the degree of dependence of δ p on the geographical variables in Central America [43]. However, the advantage of the present approach is that it uses Variance Inflation Factor to control possible multicollinearity [42].
Moran's I global statistic [44] was calculated (binary weighting, bandwidth = 4) for all three isotopic parameters, to investigate whether the residuals carry significant spatial information vital for the derivation of the desired isoscapes.

Variography
The basic function of geostatistics, the variogram, was used to explore the spatial autocorrelation structure of δ p and d-excess in the Iberian Peninsula. The empirical semivariogram might be calculated while using the Matheron algorithm [45], where γ(h) is the semivariogram and Z(x) and Z(x + h) are the values of a parameter sampled at a planar distance |h| from each other The most important properties of the semivariogram are the nugget, quantifying the variance at the sampling location (including information regarding the error of the sampling), the sill (c), that is, the level at which the variogram stabilizes, and the range (a), which is the distance within which the samples have an influence on each other and beyond which they are uncorrelated [46]. Note here, that reported ranges in the study are planar distances (d planar ) in km (EPSG: 3857), unless otherwise reported; estimated average conversion in the region: d planar × 0.755 = d geodetic . In practice, an obtained range is usually attributed to environmental processes that act on the same scale as the range of the specific variable [46], and they are known to have a physical relationship with the assessed variable. In addition, a nugget-type variogram is obtained if the semivariogram does not have a rising part and the points of the empirical semivariogram align parallel to the abscissa. In this case, the sampling frequency is insufficient for estimating the sampling range while using variography [13].
The semivariograms that were derived from the monthly GNIP data required 10 bins and a maximum lag distance of 540 km to achieve the most uniform number of station pairs in the analysis, while, in the case of those from the extended dataset, these figures were 20 bins and a maximum lag distance of 660 km. For geostatistical modeling (e.g., kriging), theoretical semivariograms have to be used to approximate the empirical ones [47]. In practice, the spherical and Gaussian models are among the most commonly used theoretical semivariograms. These were, in fact, applied in the present study by least squares model fitting [48]. While a spherical model displays a linear behavior at the origin and the sill is reached at a, a Gaussian has a repressed slope (tangent with zero slope at the origin), and its sill is a limited value that might be approximated to an infinitesimal degree, but never reached. Thus, in the case of the Gaussian model, a practical and an effective range a e = a × √ 3 have to be determined; a e is then used for further geostatistical modeling [46,49].
As the final step in the data preprocessing, semivariogram clouds were utilized to search for additional outlying values. These then had to be omitted following a detailed investigation, resulting in the final set and number of data as the input for the variography (Figure 2). The amount weighted seasonal averages were only calculated after this step has been taken; hence, monthly values that were omitted on the basis of the semivariogram clouds were also treated as missing isotopic parameters.

Isoscape Derivation with Kriging
The isoscape was derived with a three step residual kriging procedure [12,13,31,50]. As a first step, an initial grid (10 × 10 km) of stable isotope variation was derived in the region described by the multiple regression model of the supposedly driving geographic variables (LAT, LON, ELE; for details, see Sect. 2.3.2). Next, an interpolated (ordinary point kriging) residual grid was modelled while using the theoretical semivariogram (Section 2.3.3) that was fitted on to the residuals of the multivariate regression model. Note here that a methodological comparison testing various interpolation methods for the residual field of δ p in the central and eastern Mediterranean region reported that the best results were obtained while using ordinary kriging to derive the residual grid [50]. Finally, the initial and residual grids were summarized to obtain the final isoscape.
All of the computations were performed while using Golden Software Surfer 15, GS+ 10, and the raster [51] and lctools [52] packages in R [53]. The digital elevation model was retrieved from Shuttle Radar Topography Mission data [54] while using Global Mapper. Gimp 2.8 and MS Excel 365 were used for certain visualizations of the results.

Results and Discussion
The obtained significant (p < 0.05) global Moran' Is for each isotopic parameter range between 0.18 and 0.27, indicating positive spatial autocorrelation in the residuals. This encouraged the incorporation of the 'residual spatial patterns' in the derivation of the isoscapes with variography.

Variography Using GNIP Data
Variogram analysis was first conducted on monthly residuals of the δ p records of GNIP stations, and then the amount weighted seasonal averages (winter: October-March; summer: April-September) and the annual residuals (Table S2). With regard to δ 18 O p, three of 12 the semivariograms (two monthly, January, March see Figure 4A, and one seasonal, summer) were of the nugget-type, while, in the case of δ 2 H p , this was true for two of the 12 cases (two monthly: January, March). Turning to d-excess, eight out of 12 semivariograms (five monthly: October, November, January, February, March see Figure 4B, 2 seasonal: summer, winter, and the annual) were of the nugget-type.

Variography Using the Extended Database
There was a unique opportunity to extend the δp records of the GNIP stations with an additional dataset [35] covering a two year period (October 2004-June 2006; Figures 1 and 2). The station network of this regional ad hoc sampling campaign provided a seamless union with the coexisting GNIP data set ( Figure S1). The merged dataset had a greater amount of short distance data pairs than Unfortunately, a fair portion of the theoretical semivariograms proved to be unsatisfactory due to the lack of data from short distances. There is no month for which it can be said that there are both (i) stations sampling within a 50 km of each other and (ii) more than two stations sampling simultaneously in a 100 km distance. Worse still, in e.g., March 2003, there was no station pair sampling within a 100 km radius ( Figure 4A,B), a critically low degree of coverage. Therefore, the small variances over short distances were not observed, only large variance over long distances. Unfortunately, the fitting of theoretical semivariograms was not performed because of the low number and the lack of station pairs at the shortest distances, even in those cases where the semivariogram did outline a rising section ( Figure 4C,D). Besides the spatial deficiency of the GNIP network discussed above, the partial evaporation of raindrops [55] during the events with only a few millimeters of precipitation-typical in hot Mediterranean summers-can also result in local noise, further encumbering the geostatistical analysis of the δ p signal. The fact that unusual δ p characteristics were reported for S Europe [56,57] during the severe European heat waves of 2003 [58] seems to concur with this statement.

Variography Using the Extended Database
There was a unique opportunity to extend the δ p records of the GNIP stations with an additional dataset [35] covering a two year period (October 2004-June 2006; Figures 1 and 2). The station network of this regional ad hoc sampling campaign provided a seamless union with the coexisting GNIP data set ( Figure S1). The merged dataset had a greater amount of short distance data pairs than that from the GNIP stations. The Gaussian semivariograms that were fitted to the empirical variograms of δ 18 O p ( Figure 5A) and δ 2 H p ( Figure S5) obtained from the merged GNIP & NW dataset indicated a~510 km and a~420 km range, respectively. For the d-excess, a spherical semivariogram indicated a range of 440 km ( Figure 5B). If these ranges are considered to define the isotropic impact areas, the extended network provides full spatial coverage of the peninsula regarding all three parameters ( Figure 5 and Figure S5).
In the meanwhile, a shorter spatial range (below~100 km) seems to be present in the case of the primary isotopic variables, however this remains speculative and needs further investigation.
The determined spatial range (average~450 km in planar distance) indicated by the primary isotopic parameters and the range of d-excess refer to a substantial change in the spatial correlation structure of isotopes in precipitation at the scale of the zonal extension of the peninsula. The spatiotemporal variability of moisture contribution over the Iberian Peninsula can link this observation. Studies agree that the Atlantic Ocean and the Mediterranean Sea are the primary marine moisture sources, with a varying degree of spatial-and intra-annual dominance over the peninsula [26,59]. In addition, a characteristically different isotope-hydrological signature is observed between the Mediterranean (d-excess~16% ) and the Atlantic (d-excess~10% ) moisture sources (e.g., [28,60]). The recently documented division at~5 • W longitude separating the western part of the peninsula with no influence from the Mediterranean moisture in winter, from the eastern part where western Mediterranean moisture is prevailing during spring and summer seasons [27], is well in line with the above interpretation. The determined spatial range (average ~450 km in planar distance) indicated by the primary isotopic parameters and the range of d-excess refer to a substantial change in the spatial correlation structure of isotopes in precipitation at the scale of the zonal extension of the peninsula. The spatiotemporal variability of moisture contribution over the Iberian Peninsula can link this observation. Studies agree that the Atlantic Ocean and the Mediterranean Sea are the primary marine moisture sources, with a varying degree of spatial-and intra-annual dominance over the peninsula [26,59]. In addition, a characteristically different isotope-hydrological signature is observed between the Mediterranean (d-excess ~16‰) and the Atlantic (d-excess ~10‰) moisture sources (e.g. [28,60]). The recently documented division at ~5°W longitude separating the western part of the peninsula with no influence from the Mediterranean moisture in winter, from the eastern part where western Mediterranean moisture is prevailing during spring and summer seasons [27], is well in line with the above interpretation.

Isoscape of Precipitation Oxygen and Hydrogen Stable Isotopes for the Iberian Peninsula
The initial grid of δ 18 O p derived from the multiple regression model of the "Extended dataset" (Table S1) was combined with the residual grid that was obtained with the Gaussian semivariogram while using Equation (2) to derive the estimated values of δ 18 O p ( Figure 6A where A and B stand for zonal and meridional Web-Mercator coordinates (EPSG: 3857), respectively, C for elevation and D for the residual grid. With a similar procedure, the δ 2 H p isoscape was also derived (Table S1, Figure S6). Data products and estimated spatial errors of the global regionalized cluster-based water isotope prediction (RCWIP) model [16] are both freely available [61], so detailed numerical comparison between the Iberian subset of RCWIP model and the current regionally fitted δp isoscapes could be calculated by subtracting the present model values from the resampled (10 × 10 km) RCWIP ones. The δ 18 O p values that were indicated by the obtained isoscape range between −10.55 and −3.97% ( Figure 6A). The estimations error is < 0.5% at the stations, while the largest error (~0.75% ) is obtained in areas where station data were not available (as usually observed, see [14]) e.g., in the northern (Pyrenees) and SW part of the peninsula (Algarve) ( Figure 6B). A similar pattern is seen for δ 2 H p with values ranging from −76.6 to −21.1% ( Figure S6A) and estimated errors up to 6% ( Figure S6B).
The presented δ p isoscapes for the Iberian Peninsula can be compared to former products (models and maps) predicting/featuring the spatial distribution of amount-weighted long-term precipitation water stable isotopes across the Iberian Peninsula. A map of δ 18 O p distribution in peninsular Spain and the Balearic Islands was produced using composite monthly samples for the period 2000-2006 of 15 GNIP stations, belonging to the Spanish network for isotopes in precipitation (REVIP) in a polynomial model while employing latitude and elevation as predictors [29]. The REVIP-based map is not available in numerical format, so only the mapped features and some methodological differences can be evaluated. The most striking similarity is that both the current δ 18 O p isoscape ( Figure 6A) and the REVIP-based map [29] capture the orography-induced isotopic depletion for high reliefs. The polynomial model with geographical factors can be regarded as an analogue of the initial grid of this study. Thus, the added value of the presented isoscape is the consideration of the spatial variability that is retained in the model residuals irrefutably carrying sub-regional information. For example, zonal heterogeneity of vapor sources [13,14,31], which is especially interesting information from an isotope-hydrometeorological perspective and must be taken into account to improve regional isotopic landscapes. Indeed, ignoring the δ 18 O p residual variance might inspire the debatable conclusion that Atlantic Ocean is the main source of water vapor over the region [29], since the isotope-hydrometeorological dichotomy reflecting the dual moisture sources, clearly also expressed by recent atmospheric transport models [26,27], could be captured in the residual field. Finally, technical advancements in the current δ 18 O p isoscape when compared to the REVIP-based map are that the extended station network (14 more stations from Portugal, 6 from Spain and another one from Gibraltar) definitely (i) improved the accuracy of both the fitting of the regression model and the subsequent interpolation and (ii) extended the coverage of the mapping to the entirety of the peninsula.
Data products and estimated spatial errors of the global regionalized cluster-based water isotope prediction (RCWIP) model [16] are both freely available [61], so detailed numerical comparison between the Iberian subset of RCWIP model and the current regionally fitted δ p isoscapes could be calculated by subtracting the present model values from the resampled (10 × 10 km) RCWIP ones. Their difference maps (Figure 7 and Figure S7) indicate more positive values that were obtained by the current model over the mountains: such as the Pyrenees, the Cantabrian Massifs, in a small patch covering the peak region of the Beatic Cordilleras, and along the Iberian System. Major regional patterns in the difference maps are (i) lower values in the RCWIP over extended part of the Atlantic and Cantabrian regions (difference < −2.5% for δ 18 O p and < −5% for δ 2 H p ), while (ii) positive differences along the southernmost part of the peninsula (difference > 1.5% for δ 18 O p and > 10% for δ 2 H p ) indicate that RCWIP model predicted higher values south from the Beatic Range (Alboran Coast). Over a remarkably large part of the central (e.g., Meseta, Andalusia) and northeastern part (e.g., Ebro Basin, including the coast of the Balearic Sea) of the peninsula the predicted δ p values agree (Figure 7 and Figure S7). Two climatic zone domains in the RCWIP model cover the Iberian Peninsula [16], however their boundary does not separate the positive and negative field seen in the difference maps (Figure 7 and Figure S7). Moreover, the predicted values agree very well for extended areas within both RCWIP zones, so the regional scale differences likely do not reflect the RCWIP zones. The differences cannot be explained by the incorporation of the additional samples from the northwestern sector of the peninsula [35], either, because (i) the region of negative difference also extends to the Cantabrian Coast, from where no station was included in the complementary dataset and (ii) the field of positive anomalies along the Alboran Coast is situated >500 km form the area covered by the complementary dataset. It is suggested that the regional patterns seen in the difference map are evidence of a considerable improvement in the current regional prediction of the spatial distribution of long-term mean δ p across the Iberian Peninsula. northwestern sector of the peninsula [35], either, because (i) the region of negative difference also extends to the Cantabrian Coast, from where no station was included in the complementary dataset and (ii) the field of positive anomalies along the Alboran Coast is situated >500 km form the area covered by the complementary dataset. It is suggested that the regional patterns seen in the difference map are evidence of a considerable improvement in the current regional prediction of the spatial distribution of long-term mean δp across the Iberian Peninsula. The kriging error of the presented Iberian precipitation isoscape (as seen from the in-set histograms; Figure 6B and Figure S6B) is smaller by more than a factor of two for the Iberian Peninsula when compared to the RCWIP model [16], where it ranged between 0.96-0.98‰ for δ 18 Op and 7.9-8.1‰ for δ 2 Hp. Thus, the presented isoscape provides a more precise estimation of δp than the regional subset of the RCWIP model. This is in harmony with the statement, where the derivation of The kriging error of the presented Iberian precipitation isoscape (as seen from the in-set histograms; Figure 6B and Figure S6B) is smaller by more than a factor of two for the Iberian Peninsula when compared to the RCWIP model [16], where it ranged between 0.96-0.98% for δ 18 O p and 7.9-8.1% for δ 2 H p . Thus, the presented isoscape provides a more precise estimation of δ p than the regional subset of the RCWIP model. This is in harmony with the statement, where the derivation of regional isoscapes is suggested to better estimate the regional variability of δ p [16], since these do take local geographical features or meteorological parameters into account, which could influence local rainfall patterns leading to a good agreement between observed and modelled data, unlike global isoscapes [11].
It should be noted that these differences in mean predictions might partly reflect different temporal coverage. The RCWIP is a long-term climatological average that is based on pooled non-continuous decadal-scale data for 1960-2009 for δ p , while the Iberian isoscape covers a narrower interval (2004)(2005)(2006). Moreover, the differences in prediction uncertainty are almost certainly affected by different time-averaging i.e. the uncertainty in RCWIP model parameters includes uncertainty that is associated with predicting average values at stations that sample different periods of time, while, in the case of the present Iberian isoscape, this was not the case.

Conclusions
The monthly GNIP δ p records (October 2002-September 2003) for 32 stations distributed across the Iberian Peninsula appear to be unsatisfactory for use in the determination of a sampling range of isotope hydrometeorological parameters with variography. However, extended by the addition of a 21-station regional dataset (for the period October 2004-June 2006), a much denser spatial coverage was obtained, proving to be suitable for the exploration of the spatial autocorrelation structure of δ p across the Iberian Peninsula. Gaussian semivariograms were modelled for δ 2 H p and δ 18 O p , and an additional one for d-excess with ranges corresponding to~450 km in planar distance (~340 km in geodetic distance). The~450 km spatial range is thought to be related to hydroclimatological processes prevailing on the scale of the zonal extension of the peninsula and they could be connected to a spatiotemporal switch between Atlantic and Mediterranean moisture sources. A sparser network (e.g., GNIP) might also prove to be representative since the obtained spatial range of the δ p monitoring network provides a coverage for the peninsula. It means that the GNIP monitoring network is suitable for exploration of large-scale isotope hydrological processes/phenomena in the peninsula, but its verification could only be done with a denser network. Finally, the results encourage the development of precipitation stable isotope models at a sub-continental scale in further regions.