Mountain environments are extremely influenced by climate change but are also often affected by the lack of long and high-quality meteorological data, especially in glaciated areas, which limits the ability to investigate the acting processes at local scale. For this reason, we checked a method to reconstruct high-resolution spatial distribution and temporal evolution of precipitation. The study area is centred on the Forni Glacier area (Central Italian Alps), where an automatic weather station is present since 2005. We set up a model based on monthly homogenised precipitation series and we spatialised climatologies and anomalies on a 30-arc-second-resolution DEM, using Local Weighted Linear Regression (LWLR) and Regression Kriging (RK) of precipitation versus elevation, in order to test the most suitable approach for this complex terrain area. The comparison shows that LWLR has a better reconstruction ability for winter while RK slightly prevails during summer. The results of precipitation spatialisation were compared with station observations and with data collected at the weather station on Forni Glacier, which were not used to calibrate the model. A very good agreement between observed and modelled precipitation records was pointed out for most station sites. The agreement is lower, but encouraging, for Forni Glacier station data.

1. Introduction

Responses of the mountain environment to climate change represent one of the most studied topics in recent years [13]. The evident and fast landscape modifications occurring in the last decades testify the vulnerability of mountain areas to climate change (e.g., [4]). At higher altitudes cryosphere is shrinking and in the last 60 years Italian Alpine glaciers have lost about 30% of their volume [5], with a huge reduction of glacier surfaces [6] and a progressive widening of proglacial areas, which lead to new colonisation of glacier forelands [7, 8] and debris covered surfaces by vegetation [911] and animals [12, 13]. Such geomorphological changes involve not only the abiotic and biotic components but also the landscape fruition as they modify the local geoheritage [14], geodiversity [15], hazard, and risk scenarios [1618].

As glacial and geomorphological processes are strictly climate related, accurate reconstruction and analysis of present and past climatic conditions are crucial. Qualitative reconstruction of the past climatic characteristics over long time scales is possible thanks to geomorphological and biological paleoclimatic indicators (e.g., typical features of glacial and periglacial environment, debris covered glaciers and rock glaciers, tree remnants under glacial deposits, and pollens) [1922]. Quantitative reconstructions, instead, come from dendroclimatic analysis [2325] or from meteorological observations that in Italy began to be collected regularly in the last decades of the XVIII century [26]. The longest records of meteorological data are however generally available for anthropized areas (in Italy they concern Milan, Padua, and Turin), with some excellences in high mountain environment, such as Capanna Margherita (Punta Gnifetti 4554 m, since 1899), Sonnblick (Austria, 3106 m, since 1886), and Jungfraujoch (Switzerland, 3466 m, since 1930). In the last decades many automatic weather stations (AWS) have been positioned over glaciers and other high-elevated sites (e.g., at Forni Glacier in Valtellina (SO) [27]).

Therefore, mountain environments, especially glacier forelands, are extremely influenced by climate change but are also often affected by the lack of long and high-quality meteorological data, which limit the ability to deeply investigate the acting processes at local scale and their influence/impact on the biotic and abiotic environments as well as on the human component. This problem may be reduced, at least partially, by methodologies allowing us to extrapolate the existing observations over such areas. Many interpolation methods have been developed so far and allow reconstructing the high-resolution spatial distribution and temporal evolution of a meteorological variable over an extended area by integrating punctual station data and the topographic features of the surface, extracted from Digital Elevation Models (DEM). These methodologies allow us to obtain the climatic signal for a number of points much higher than the sparse station network and to get information even for those limited areas where no meteorological stations are present [28]. Results can be moreover integrated with dendroclimatic reconstructions, which can go back for hundreds to thousands years [29]. However, the ability of the model to reconstruct past climate series strongly depends on data distribution and availability as well as on orographic complexity of the area under investigation.

Temperature, precipitation, solar radiation, and wind are the main parameters that need to be widespread on glacier forelands. The present work is focused on precipitation as its distribution, characteristics, and trend have a great influence as triggering factors for hazard [3, 30]. They can represent limiting factors for biological community [31] and they can affect human fruition of the higher areas [32]. Moreover, the ability of quantifying the precipitation distribution (together with temperature) and its features over Alpine regions is extremely important to improve the estimation of mass balances over glaciers or hydrological balances over catchments.

Within this context, the paper aims at (I) presenting a methodology to reconstruct monthly precipitation time series (1913–2015) for any point of a 30-arc-second-resolution grid, covering a complex mountain area centred on the Forni Glacier (Central Alps) where an AWS is located since 2005; (II) analysing the spatial distribution of precipitation normals and precipitation temporal trends; (III) comparing the modelled precipitation series over Forni Glacier with the AWS records; and (IV) discussing the possible applications to multidisciplinary researchers focusing on glacier foreland changes.

2. Materials and Methods

2.1. Study Area and Data

The area selected for the present work is centred on the Forni Glacier, one of the largest Italian valley glaciers (11.34 km2 [5]), and it covers about 2000 km2 (grey shaded area in Figure 1).

In the last decades, this mountain area, as well as the similar ones all over the Alps, has suffered a rapid and deep modification due to climate change [33]. The most evident result is a generalized shrinkage of glacier bodies and, consequently, the expansion of glacier forelands with changes in the morphology of the upper valley portions. A deeply known moraine system (see for details [21]) characterises the Forni Valley allowing us to precisely date the glacier advancing phases since the Little Ice Age [34, 35] and to reconstruct the progressive glacier retreat and the following tree colonisation (ecesis) [8]. The germination of new plants in the proglacial areas, together with the increase of the glacial debris coverage, due mainly to gravity processes and ice ablation, contributes significantly to the deep transformation of the local glacial landscape. The debris coverage over the Forni Glacier ablation tongue is estimated to be increased from 26.7% in 2003 to 48.1% in 2015 [36]. This variation is in turn responsible for local microclimatic changes also in relation to the type of lithology constituting the debris. Land cover change caused by glacier shrinkage can modify climate through different mechanisms, some directly perturbing the Earth radiation budget and some perturbing other processes [37]. In fact, land cover change can alter the surface energy and moisture budgets through changes in evaporation and the fluxes of latent and sensible heat, directly affecting precipitation and atmospheric circulation as well as temperature [37].

The study area includes the Upper Valtellina (11 stations), the western part of the Venosta Valley (9 stations), and the western part of Sole Valley (9 stations), located in the Central Italian Alps, and a portion of the Grisons Canton (1 station), on the Swiss side of the Alps. A great part of this region is included within Stelvio National Park; see Figure 5 and Table 4. This area is characterised by a significant orographic heterogeneity with elevation values ranging from approximately 1000 m a.s.l. of the valley bottom to the 3769 m a.s.l. of the highest peak (Mount Cevedale). A bigger area was also considered (black border square region in Figure 1) in order to strength model computation of precipitation fields on the inner area. The considered weather stations (grey shaded area in Figure 1) can be grouped according to their distribution with respect to the relief features. Two stations belong to the high mountain environment, whose lower limit roughly corresponds to the upper limit of vegetation [38]; 8 stations are located on slopes at middle altitude (below the treeline) and 6 more stations are located at Passes and plateau. The remaining 14 stations are positioned at the valley bottoms (Table 4).

The database considered in the present work is composed by 30 daily precipitation series, covering the study area, which were retrieved from several Italian regional and subregional sources (Meteotrentino for the autonomous province of Trento, the Hydrological Office of the autonomous province of Bolzano/Bozen, the weather service of Regional Agency for Environmental Protection ARPA, in Lombardia, and the Geological Monitoring Service, CMG, weather network also in Lombardia) and from MeteoSwiss. Details on data sources are listed in Table 1. The series cover different time intervals in the 1913–2015 period and their elevation ranges are not equally distributed (Figures 2(a) and 2(b)). Fourteen of those series have been obtained by merging records from the traditional networks of mechanical stations with records from the new networks of automatic stations. This merging was performed only if distance between the old and the new site was within 100 m and elevation difference within 30 m. Moreover, in some cases the two series had an overlapping period that supported their homogeneity assessment.

The daily data were subjected to a preliminary gross-check. Then, they were subjected to a gap-filling procedure in order to minimise the effect of sparse gaps in the daily records on the computation of monthly totals. For each station (test), all months with at least 15 available daily data were subjected to daily fulfilment procedure, where each daily reconstruction was performed by considering the daily records of a reference station. The reference station was selected considering distance, elevation gap, common days in the considered month, and correlation value with the test series. Correlation was computed by Pearson method and all series with correlation values less than 0.7 were discarded. Therefore, reference station was the nearest one showing the highest correlation and number of common data. The missing daily values were reconstructed using the multiplicative anomaly method [39]. More in detail, the average daily precipitation was computed for test and for reference stations on the basis of their common period, with these two averages, it was possible to compute the ratio that worked as scaling factor. The missing day in test series was reconstructed by rescaling the corresponding day in reference station using the aforementioned scaling factor.

On average, the percentage of daily data reconstructed for each station was about 1% (the maximum was about 5% at station Arnoga Valdidentro). The robustness of the filling procedure was evaluated by applying it to the available known records for each daily series and comparing the measured values to the simulated ones in terms of both daily precipitation values and number of wet days (Figure 3). Uncertainties were evaluated for the reconstruction of both absolute daily precipitation and number of wet days (days with precipitation greater than or equal to 1 mm) for each month. The uncertainties showed that for absolute daily cumulate precipitation root mean square error (RMSE) was around 2.5 mm and Mean Absolute Error (MAE) was 1.0 mm, with a total of 449534 simulated days. MAE was found to be more affected by the nonrainy days (precipitation lower than 1 mm), so we decomposed the computation of those errors into several rain classes. The chosen classes and the corresponding errors are summarised in Table 2 and highlight the increase in MAE and RMSE values with increasing precipitation, whereas the reconstruction of monthly wet days provided a quite good result: on a total of 14912 simulated months MAE was 1.5 days and RMSE was 2.1 days. These results are also reported as scatter plots in Figures 3(a) and 3(b).

After the daily fulfilment, monthly dataset was computed. Whenever missing daily data were still present, the corresponding monthly cumulate was not evaluated. The increment of daily data, even though it was little, had permitted to increase significantly the monthly precipitation availability: on average the benefit was around 15%, and the maximum increment was obtained for station Arnoga Valdidentro, that passed from 77 to 182 available months.

In order to investigate monthly series in terms of quality and homogeneity and to deal with the points at the borders of study area in the climatological reconstruction, they were integrated by several monthly precipitation series retrieved from a quality-checked database recently set up to compute the 1961–1990 monthly precipitation climatologies for Italy [40]. All the available stations in a wider area (latitude north 45°00′–47°48′ and longitude east 8°36′–12°00′) were included reaching a total of 734 series and a spatial density of about 9 stations/103 km2.

Besides the 30 station series, we also considered, but did not include in the database used as input for the model, the daily observations collected by an AWS (named AWS1 Forni) that has been operating on the ablation tongue of the Forni Glacier since 2005 and monitors continuously the most important meteorological variables. It represents a precious and unique source of data directly measured on a glacier surface and they were used as validation of the model outputs.

2.2. Monthly Quality-Check and Homogenisation

The quality-check and homogenisation activities on the 30 monthly series in the domain were carried out in order to remove outliers and spurious values and to detect nonclimatic signals due to, for example, changes in station location, station malfunctions, changes in surroundings, sensors, and/or measurements techniques [41].

Quality-check was performed for all the 30 stations inside the domain, by following the method described in Crespi et al. [40]. It is based on the comparison between the measured series (test) and the simulated one reconstructed on the basis of ten neighbouring stations (reference) with enough data in common with the station under consideration.

The homogenisation algorithm is based on the Craddock test [42] and the Expert Judgment method [43, 44]. It is based on the assumption that if two series are homogeneous their ratios should be constant in time. Therefore, for each test series, five reference stations are chosen based on distance, altitude criterion, and amount of common data. The monthly series of cumulative relative differences (Craddock series) between test and each reference were computed, and inhomogeneous periods were selected manually and corrected by the proper correcting factors; nine series out of 30 were homogenised for relative short periods. In Table 3, the homogenised series are listed together with the corrected periods, while in Figure 4 the annual series of Bormio station before and after the homogenisation is reported as an example.

After these procedures, the monthly normals in the period 1961–1990 were computed for each station. In order not to introduce biases in the normal computation due to unequal number of available data, missing months were reconstructed by following the procedure described in Crespi et al. [40] and based on the same method applied for quality-check.

2.3. Interpolation Models and Anomaly Method

The 1961–1990 high-resolution precipitation climatologies were computed for each grid cell of a 30-arc-second-resolution DEM (GTOPO 30 DEM, U.S. Geological Survey) covering the whole study domain. This DEM has a root mean square error (RMSE) of about 18 m [45]. Since many climatological works (see, e.g., [46]) proved that precipitation distribution is strongly linked to orography, especially elevation, and the interpolation techniques taking into account terrain features are proved to provide better performances, the monthly normals were distributed applying a Local Weighted Linear Regression (LWLR) of precipitation versus elevation [47]. LWLR computed the precipitation normal at each DEM cell by considering its elevation :where and are the coefficients of the weighted linear precipitation-elevation regression performed at the considered cell by selecting the stations with the highest weights.

Station weights () are expressed as the product of several weighting factors in the form of Gaussian functions (2) in which the distances and the level of similarity between the station cells and the considered grid cell in terms of orographic features (i.e., elevation, slope steepness, and slope orientation) are taken into account:where is the difference of the parameter values between the station and the grid cell and is the coefficient regulating weight decrease.

Orographic features of each station were extracted from the closest DEM cell, while the most suitable decreasing coefficients to be used for the weighting factors were optimised month-by-month following the procedure described in Crespi et al. [40].

In order to assess the suitability of the chosen interpolation model, we compared it with the widely used Regression Kriging (RK) approach. This comparison allowed us to test model robustness and to assess the most suitable spatial scale at which the precipitation-elevation relationship has to be studied on the domain.

In RK, a global precipitation-elevation regression is computed on all the stations in the domain and the station residuals are then interpolated onto high-resolution grid by ordinary kriging. Precipitation normal at cell is finally obtained as follows:where and are the global regression coefficients, are the kriging weights, and are station residuals. Further details on RK approach can be found in, for example, Hengl et al. [48].

Moreover, in order to define the actual spatial scale at which the direct effects of elevation on precipitation are important [49], a Gaussian filter was applied to smooth out terrain features, while retaining the 30-arc-second-resolution. The smoothing of the DEM was performed by assigning to each cell an elevation obtained as a weighted average of the elevations of the surrounding cells, with weights provided by a Gaussian function decreasing to 0.5 at a distance of 3 grid steps from the cell itself; the degree of smoothing was chosen after an optimisation procedure (see more details in Crespi et al. [40]). The DEM is depicted in Figure 5. Climatologies were therefore computed on this smoothed version of the DEM (3smDEM), which is also used to assign to the stations the orographic parameters required in the interpolation procedure.

After computing the distribution of 1961–1990 precipitation normals on the domain, the century records for the period 1913–2015 were obtained for each grid point by applying the anomaly method [50]. This methodology is based on the assumption that the spatiotemporal behaviour of meteorological variables can be described by the superimposition of a constant field (i.e., the climatologies) and the departures from them (i.e., the anomalies). At this aim, long and high-quality homogenised series are required to avoid reconstructing nonclimatic signals. Due to the higher spatial coherence of anomalies in respect to climatologies, interpolation approaches considering larger scales and neglecting local surface features are preferable [28]. Monthly precipitation series from all stations were converted into anomaly series as the ratio to the corresponding 1961–1990 monthly climatologies (normals). Afterward they were interpolated onto the 30-arc sec grid by averaging the anomalies of neighbouring stations (within 200 km from the grid point) by means of Gaussian weighting functions (in the form of (2)), taking into account their distance from the considered grid cell and elevation difference.

Due to the evolution of data coverage during the considered time interval, the weight decrease was regulated for each year according to station availability. While for elevation weight a constant halving coefficient was set, for distance weight it coincided with the mean radius from the grid point within at least 3 available stations that were included. This value ranged between 45 km at the beginning of the period to 10 km for the central decades when station coverage was maximum. High-resolution field of 1913–2015 monthly precipitation series in absolute values was finally obtained by multiplying gridded anomalies for 30-arc-second-resolution climatologies.

2.4. AWS1 Forni as Model Validator

The statistical model, as previously described, was set up using a net of weather stations ranging from valley bottom to the mountain top. Each station insists on a specific area and altitude gap, and its influence decreases progressively (see Section 2.3). Each station represents its surroundings, especially in such a complex terrain, and overlaps to the predominant climate signal a microclimate forcing. In order to assess the ability of the reconstruction approach to extrapolate the local microclimatic conditions of the area, especially over those points which are less monitored by in situ measurements, we compared the modelled series for a gird point over the Forni Glacier with the precipitation observations collected by AWS1 Forni. This station was not included in the database entering in the model, so its data were completely new and a very good start point for the validation.

The station is equipped with a not-heated rain gauge, which measures the summer (from June to October) liquid precipitation, and two sonic ranger sensors, which measure the snow depth. For the estimation of daily snow water equivalent (SWE) an approach based on hourly snow depth data and a fresh snow density of 140 kg/m3 was applied [27, 51]. This method was validated in a previous study by comparing the estimated SWE values with data measured by means of manual snow pits and recorded by an automatic snow pillow [52].

Until October 2014, the AWS1 Forni was situated on the lower sector of the eastern tongue of the Forni Glacier. Because of the increased number of crevasses that made it very difficult and hazardous to reach the station site, the AWS was moved to an upper area. Unfortunately, due to the formation of ring faults that could compromise the stability of the station itself [53, 54], in November 2015 the AWS was moved on the central tongue. All the changes in AWS location are reported in Table 5. The actual AWS1 Forni position was set to the average of the three locations. The nearest DEM cell to this position is 2101 (Figure 6), which was used together with the other eight nearest points to compare the model with the station data.

3. Results and Discussion

3.1. Climatologies

The monthly climatologies (1961–1990 normals) were evaluated by LWLR for the inner domain reported in Figure 1 and model accuracy was evaluated in terms of Leave-One-Out (LOO) reconstruction of monthly normals of the 30 stations of this area, that is, by excluding the station meant to be reconstructed in order to avoid self-influence. Model uncertainties were estimated in terms of monthly MAE, RMSE, and BIAS. In order to assess the suitability of the chosen approach for the peculiar study area, model results were compared to those obtained by RK and the corresponding monthly errors are reported in Table 6. On average LWLR had shown lower values for MAE and RMSE with respect to RK, except for summer months when RK performances were comparable or even better than LWLR. The reason could be found in the peculiar LWLR extrapolation ability: the local precipitation-elevation regression tends to enhance precipitation amounts in summer when the highest station monthly normals occur. However, the annual evolution of climatologies for each considered station was well captured by both models.

Regarding the annual climatological map depicted in Figure 7, a very dry zone was reconstructed in Venosta Valley (northeastern part of the box), with less than 400 mm/yr, while the mountain chain of Ortles-Cevedale is more wet, with 1000–1100 mm/yr. Meanwhile the Upper Valtellina (Bormio, Valfurva, and Livigno) has an intermediate precipitation pattern (around 750–900 mm/yr). Moreover, in the western part of the area is visible a very high precipitation region (up to 1400–1500 mm/yr): this is related to the high precipitation measured across the Bernina Pass that has Stau conditions from both north and south sides. Considering an average on the entire area the rainfall is 880 mm/yr.

As pointed out by annual precipitation distribution, Ortles-Cevedale chain acts as a pluviometric shadow, as it stops the southerly wet air masses. This behaviour is even more evident observing the seasonal climatologies in Figure 8, especially during winter. Moreover the annual cycle points out the contrast between a very dry winter (even less than 100 mm/season, especially on Upper Valtellina and Venosta Valley) and a very wet summer (up to 400 mm/season, especially on Ortles-Cevedale chain). The relevant seasonality of the shadowing effect of the Ortles-Cevedale chain is a consequence of the fact that atmospheric circulation during rainy days exhibits strong differences along the year. Considering an average on the entire grid the amount of precipitation is 124 mm for the winter season, 228 mm for spring, 299 mm for summer, and 230 mm for autumn with minimum in February (40 mm) and maximum in August (107 mm). The lower winter values can be explained by the typical atmospheric circulation which is present in winter when precipitation occurs over North-Eastern Italy: it is characterised by southeasterly advection blowing wet air from the Adriatic Sea. Those wet air masses, blocked by the Central Alps, discharge their moisture over the South-Eastern side of Ortles-Cevedale chain. In the meanwhile, the North-Eastern part of the domain features analogous conditions due to the presence of Swiss Alps, which block frontal systems moving from the north and force them to lose humidity on the Switzerland side. The summer season, on a meteorological point of view, is characterised by a marked contrast between warm and moist air coming from the Po Plain and a relatively fresh continental air in the north side. Within this configuration, convective cells are promoted and convective precipitation is developed, in contrast with winter season where precipitation is generated mainly by orographic barring.

The same maps of Figures 7 and 8 were also built by the RK (not shown). While the annual cycle and the overall spatial distribution are comparable to those reconstructed by LWLR, pluviometric gradients were more smoothed out due to the global regression approach applied. LWLR local approach allows us to obtain a very detailed spatial distribution of the interpolated variable, and this could be particularly interesting for heterogeneous orographic regions, such as the considered one, but may lead to artefacts in case of poor station coverage over complex terrains. On the contrary, by considering elevation-precipitation relationship at a larger scale RK smooths more out the punctual climatological signal but provides more robust results in cases of low data availability.

3.2. Anomalies and Trends

After reconstructing the climatological field, the 1913–2015 gridded anomaly dataset was computed, and its suitability for the analysis of the time variability in precipitation series for a single grid point of interest or for the whole domain was assessed. Model ability in estimating anomaly series was evaluated in terms of BIAS, MAE, and RMSE values obtained by the reconstruction of station anomalies in LOO approach. Monthly MAE ranged between 0.28 for winter months and 0.14 for summer months.

The evolution of MAE during the entire period was also analysed in order to evaluate the reconstruction reliability. Due to the changes in quality and distribution of stations during the period, we obtained higher values of MAE, around 0.30, at the beginning and at the end of period, related to both lower data coverage and quality. In the middle part of our dataset MAE values were about 0.20. The evolution of BIAS was more uniform, especially since 1960, and approximately around zero.

The rather low LOO-errors of the simulated anomalies are due to rather high common variances between the anomaly records of the different stations: the scatter plot of correlation versus distance (Figure 11(b)) gives evidence of correlation coefficients over 0.7 even for the highest possible distances in the study area. The high common variance between the station anomaly records is also highlighted by Principal Component Analysis (PCA), the first eigenvector explaining more than 80% of the variance of the data set of inner area stations.

The high spatial coherence of the station anomaly records suggests investigating trends focusing only on an average inner area anomaly record. Short- and long-term trends for the whole domain were evaluated both on annual and seasonal scale. Long-term trends were computed using Theil-Sen test, and their statistical significance was proved by Mann-Kendall estimator [5557]. Even if slope was slightly negative for winter and spring series, suggesting a centennial reduction of precipitation over the area, and positive for summer, autumn, and annual scale, the long-term trend has shown a statistical significance below the 95% confidence level. Short-term trends could be highlighted using an 11-year Gaussian window filter with a 3-year standard deviation (Figure 9) which smooths out random variability and evidences more clearly the decadal fluctuations. In order to better describe the complex time evolution of climatic signal, which cannot be fully captured by the Theil-Sen slope on 1913–2015 period, a running-trend analysis was also performed on annual and monthly scale. Theil-Sen slopes and significance levels were estimated within windows whose widths range from 20 years up to the entire series length, running from the beginning to the end of the series [26]. Figure 10 shows the results of the running trend performed on annual precipitation anomalies, where colours represent the values of trend while size of pixels describes its significance. As already pointed out by long-term trend analysis, trend values and significance are mostly negligible over the longest time intervals, while more evident trends occurred considering shorter time windows. Years between 1925 and 1945 featured a clear negative trend, while annual precipitation amount was found to increase in both 1940–1960 and 1960–1980 periods.

3.3. Cross-Validation and Comparison at AWS1 Forni

The inner area stations reconstructed precipitation records—obtained superimposing the climatologies and the anomaly records—were also validated by means of the LOO approach. The MAE turned out to be about 15 mm/month. The lowest value was obtained for Silandro (ID 16 in Table 4) with less than 9 mm/month, while the highest value was found for Valdisotto (ID 26 in Table 4) with over 30 mm/month. The MAE errors of all inner area stations are represented in Figure 5.

We also computed the linear regression between the simulated and measured values at all inner area station sites (Figure 11(a)) and we foundThe highest common variance between the reconstructed and the observed record was found for Rabbi Somrabbi (ID 9 in Table 4); the lowest was found for Valdisotto.

A seasonal analysis was also performed using two macro-seasons: fall-winter (from September to February) and spring-summer (from March to August). The results were similar to the all-data regression. They are reported below:Spring-summer reconstructions were generally less accurate than fall-winter ones and this could be related to the lower spatial coherence of precipitation in the summer period.

Finally, the model performances were further assessed by comparing modelled rainfall series and the measured values collected by the AWS on Forni Glacier, not included in the input dataset and located in a very remote and peculiar area. The precipitation series from the AWS on Forni Glacier counts 3563 days, from 1 October 2005 to 30 September 2016, and the number of complete months is 97 (no missing days). All trustworthy monthly data between 2005 and 2015 were compared with the modelled series of grid point 2101, which is the nearest to the average position of the AWS, and the computed errors are 33 mm/month for the MAE and 42 mm/month for the RMSE. The comparison of monthly data is shown in Figure 12(a). Considering an average of cell 2101 and its eight nearest neighbours, the errors raise to 38 mm/month of MAE and 48 mm/month of RMSE; this could be due to the inclusion of lower values far from the peak precipitation area of the watershed. These rather high errors depend indeed on the difficulties in capturing the exact position of the transition between the wet areas influenced by southeasterly currents and the very dry Forni Valley; however, they could also be strongly influenced by the high uncertainties of the SWE evaluation method. In fact, comparing the SWE values with the ones measured for two winter seasons by a snow pillow, a RMSE value of 45 mm w.e. was found [52]. In addition, even if the snow surface over the Forni Glacier seemed to be homogenous, a large spatial variability of snow depth was observed carrying out numerous measurements through a snow weighting tube (Enel-Valtecne©) around the AWS in February 2015 with a standard deviation of 29 cm of snow; see, for example, Senese et al. [52]. This is due to the high roughness of the ice surface that causes different snow accumulation rate over the glacier surface. Consequently, the error between the modelled series of grid point 2101 and the ones observed at the AWS1 Forni could be enhanced even by the high spatial variability of the snow depth in high mountain areas.

A deeper temporal analysis highlighted a high seasonality between simulated and measured precipitation (Figure 13); in fact the winter period is underestimated by the model; meanwhile the summer period is quite well depicted. The uncertainties follow the same pattern: December to February RMSE (DJF, meteorological winter) is 59 mm/month, and MAE is 54 mm/month, instead June to August RMSE (JJA, meteorological summer) is 24 mm/month, and MAE is 18 mm/month. The divergence between simulated and measured data could be explained by the peculiar precipitation distribution over the area. The orographic effect of Ortles-Cevedale chain is captured in the reconstructed field (Figure 8), but discrepancies with AWS1 Forni data suggest that the actual position of regime transition could not be fully resolved by the model together with the chosen DEM grid spacing. In fact, LWLR places the maximum rainfall (snowfall) over the crest, but especially for snowfalls, the wind could transport snowflakes for several kilometres, and it could create deposit far away from the precipitation maximum. Other possible sources of errors could be the conversion of snow height into water equivalent (see Section 2.4), for which a constant snow density was used, and the measurement itself of snow depth. In fact, these measurements could be overestimated due to the wind that acts on snow cover, changing its distribution and accumulating snow under the sonic ranger.

In order to investigate wind effect, the wind patterns were compared to days with big discrepancies between simulated and measured data over all the period. Wind records were retrieved from the anemometer installed at AWS1 Forni.

Daily differences were calculated as BIAS between observed and simulated rain, and minimum threshold of 20 mm was chosen. However, as shown in Figure 12(b), many points lay under 5 m/s of wind speed, and only eight underestimated days have a mean wind speed greater than 5 m/s. Only the 16% of points exceeded 5 m/s threshold, and they provided no evidence of a significant relationship occurring between precipitation underestimation and wind speed.

4. Conclusions

A statistical approach to project precipitation data on a high-resolution grid and based on anomaly method was applied over a complex mountainous terrain in the Central Italian Alps. The database counts up to 734 rain gauges, and the area where the climatologies were computed includes 30 quality-checked and homogenised stations. The 1961–1990 precipitation climatologies were spatialised over the study domain by LWLR starting from the assumption that a strong relationship occurs at local level between precipitation and orographic features. LOO model errors ranged between 4 mm and 10 mm, in terms of mean monthly MAE. These performances were compared to those obtained by global RK that had a worse winter and better summer performances, with a MAE ranging between 6 mm and 11 mm. Those results allowed us to assess the goodness of the proposed method for the reconstruction of monthly precipitation time series at a specific point of the high-resolution grid. We also appreciated the different spatialisation properties of the two approaches. In fact, LWLR can better represent orographic features than RK that, instead, smooths out such details; on the contrary, the first one is more susceptible to outliers and/or overestimations (especially during summer) than the second one, in particular over those areas less covered by stations.

Regarding the reconstruction ability of 1913–2015 anomalies, monthly MAE along the whole period ranged from 0.14 to 0.28.

We verified the absence of statistically significant long-term trends for the considered area; however short-term trends were present. In particular, three positive tendencies were found (1940–1960, 1960–1980, and the last decade) and one negative during the 1920s. The spatialised climatologies revealed a very dry area in the Northern part of our study domain, corresponding to Venosta Valley; by observing seasonal climatologies, a two-pronged behaviour of precipitation for Upper Valtellina was found: a very dry winter, where mountain barring mostly acts, opposed to a very wet summer, when convective precipitation is prevalent.

The cross-validation performed for the 30 inner area stations showed that the key factor for reconstructing unbiased precipitation fields is correctly capturing the spatial variability of the monthly precipitation normals. It is therefore very important to exploit the existing data as much as possible in order to have the best spatial resolution of these data and to recover all possible records even if they cover only rather short periods. On the contrary, the high spatial coherence of the anomalies allows capturing the temporal evolution of precipitation over the investigated area with a much lower number of stations. In this case, however, homogeneity is a key issue and the availability of a reasonable number of records turns out to be very important for homogenisation procedures.

The spatialised fields were also compared with the in situ measurements of the supraglacial automatic weather station located since 2005 at Forni Glacier, which was not included into the model database. We have obtained that the model had well depicted the main pluviometric features of the area, but it faced greater difficulties during winter, when the Upper Valtellina is drier with respect to Sole Valley in Trentino. The observed uncertainties could be explained with an underestimation of the model in this watershed zone for a misplacing of the transition zone, or with an overestimation of winter snow measurements and/or distribution.

In conclusion the applied approach reveals many advantageous aspects; in fact, from 30 in situ measurements it allowed reconstructing at high resolution the precipitation signal over an area of about 2000 km2. Therefore, this approach could be particularly helpful if applied to a rather complex mountainous terrain, where direct measurements are poor and quite short, in order to understand and study fragile environments in the frame of climate change. However further researches are necessary to better investigate the role of precipitation, temperature, and wind on the biological component of glacier surfaces and forelands, especially if we consider that they represent a changing habitat for life.

Conflicts of Interest

The authors declare that they have no conflicts of interest.


The authors wish to thank the Regional Agency for Environment Protection of Lombardy (data series is available on the web at http://www.arpalombardia.it/siti/arpalombardia/meteo/richiesta-dati-misurati/Pagine/RichiestaDatiMisurati.aspx), the Geological Monitoring Service of Lombardy, the Weather Service of Trento Province (https://www.meteotrentino.it/dati-meteo/info-dati.aspx?ID=3), the Weather Service of Bolzano Province (http://dati.retecivica.bz.it/it/dataset/misure-meteo-e-idrografiche), and the Federal Weather Service of Switzerland (https://gate.meteoswiss.ch/idaweb/login.do) for supplying weather data. The AWS1 Forni has been developed by UNIMI, is hosted in the Stelvio Park-ERSAF, and is presently managed by UNIMI DESP with the financial support of Sanpellegrino SpA (Levissima). The acquired data are not only used for the scientific research of the UNIMI staff but also shared with the WMO (World Meteorological Organization) in the framework of some international programs (i.e., SPICE and CRYONET).