# High-Resolution Monthly Precipitation Fields (1913–2015) over a Complex Mountain Area Centred on the Forni Valley (Central Italian Alps)

Authors: Golzio Alessio; Crespi Alice; Irene Maria Bollati; Senese Antonella; Guglielmina Adele Diolaiuti; Pelfini Manuela; Maugeri Maurizio
Journal: Advances in Meteorology
 (2018)

---

## Abstract

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.


---

## Body

## 1. Introduction

Responses of the mountain environment to climate change represent one of the most studied topics in recent years [1–3]. 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 [9–11] 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 [16–18].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) [19–22]. Quantitative reconstructions, instead, come from dendroclimatic analysis [23–25] 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).Figure 1
Study area domains: grey shaded the inner box area (longitude E 10°10.98′–10°49.14′ and latitude N 46°15.6′–46°37.34′) centred on the Forni Glacier, with the solid black square bordering the outer area (longitude E 8°36′–12°00′ and latitude N 45°00′–47°48′). Black dots indicate the weather stations considered.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 Figure5 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 thehigh mountain environment, whose lower limit roughly corresponds to the upper limit of vegetation [38]; 8 stations are located on slopesat middle altitude(below the treeline) and 6 more stations are located atPasses and plateau. The remaining 14 stations are positioned at thevalley 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 Table1. 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.Table 1
Meteorological networks considered in this work. References for datasets are available at Acknowledgments section.
Network State Region, province ARPA Lombardia, Meteo Italy Lombardia, Sondrio ARPA Lombardia, CMG Italy Lombardia, Sondrio Meteotrentino Italy Trentino Provincia autonoma di Bolzano, Bozen Italy Alto Adige Meteo Svizzera, Meteo Swiss Switzerland GrigioniFigure 2
Availability of weather stations (a) from 1913 to 2015 and station distribution according to the altitude (b), with classes of 100 m. The inner study area is represented in dark green, and the outer area is represented in salmon.
(a) (b)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 (Figure3). 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).Table 2
Reconstruction of absolute precipitation skill for precipitation classes subdivision.
Rain class (mm) Days MAE (mm) RMSE (mm) [0.0–1.0) 326043 0.22 0.96 [1.0–2.5) 32270 1.42 2.10 [2.5–5.0) 27178 2.13 2.86 [5.0–7.5) 17492 2.83 3.69 [7.5–10.0) 10898 3.45 4.54 [10.0–15.0) 14115 3.97 5.26 [15.0–20.0) 8056 4.67 6.22 [20.0–30.0) 7776 5.49 7.40 [30.0–40.0) 3180 6.55 8.87 [40.0–50.0) 1344 7.56 10.31 [50.0–75.0) 984 9.16 12.29 [75.0–100.0) 159 11.92 16.24 [100.0–125.0) 30 13.98 17.48 [125.0–150.0) 7 28.93 35.38 [150.0–200.0) 2 21.67 22.11Figure 3
Scatter plot of simulated versus measured monthly wet days (a) and daily rainfall (b). The skill score is the red line (1 : 1), while blue lines represent the ratios 2 : 1 and 1 : 2.
(a) (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.Table 3
List of homogenised precipitation series and their homogenised periods.
Station name Network Homogenised periods Bormio ARPA Lombardia 1930–1959 Cancano (Valdidentro) ARPA Lombardia 1951–1969 Solda Provincia di Bolzano 1930–1965, 1966–1995 Careser (dam) Meteotrentino 1975–1992 Cogolo Pont (Centrale) Meteotrentino 1933–1940, 1990–1999 Passo del Tonale Meteotrentino 1923–1930, 1959–1971 Mezzana Meteotrentino 1990–2000, 2009–2016 Malga Mare (Centrale) Meteotrentino 1974–2016 Peio Meteotrentino 1930–1950, 1944, 1960–1972Table 4
Inner area station characteristics: the ID column refers to Figure5.
ID Station name Province Lat N Lon E Elevation [m] Positioning Side 1 Grosotto Sondrio 46.286 10.261 613 Valley bottom 2 Premadio centrale Sondrio 46.485 10.351 1305 Valley bottom 3 Trafoi Bolzano 46.550 10.508 1570 Mid mountain E 4 Pezzo Brescia 46.283 10.519 1557 Valley bottom 5 Agumes Prato Bolzano 46.621 10.582 907 Valley bottom 6 Gioveretto dam Bolzano 46.509 10.725 1851 Mid mountain SE 7 Ganda Bolzano 46.550 10.786 1257 Mid mountain NW 8 Lago Verde Bolzano 46.467 10.803 1488 Mid mountain E 9 Rabbi Somrabbi Trento 46.410 10.804 1352 Mid mountain S 10 Santa Caterina Valfurva Sondrio 46.415 10.498 1760 Valley bottom 11 Forni Sondrio 46.421 10.564 2118 Mid mountain Bottom 12 Cancano Sondrio 46.514 10.323 1855 Mid mountain Bottom 13 Bormio Sondrio 46.466 10.371 1225 Valley bottom 14 Solda Bolzano 46.514 10.594 1905 Valley bottom 15 Lasa Bolzano 46.615 10.693 865 Valley bottom 16 Silandro Bolzano 46.622 10.779 698 Valley bottom 17 San Giovanni Martello Bolzano 46.526 10.736 1616 Mid mountain E 18 Pian Palù (dam) Trento 46.337 10.614 1800 Mid mountain NE 19 Careser (dam) Trento 46.423 10.699 2600 High mountain NW 20 Mezzana Trento 46.313 10.796 905 Valley bottom W-E 21 Passo del Tonale Trento 46.263 10.597 1875 Mid mountain Pass 22 Peio Trento 46.364 10.678 1585 Mid mountain SSE 23 Sondalo, Le Prese Sondrio 46.350 10.355 950 Valley bottom 24 Aquilone Sondrio 46.401 10.355 1082 Valley bottom 25 Grosio, Fusino dam Sondrio 46.327 10.246 1220 Mid mountain S 26 Valdisotto Sondrio 46.384 10.330 2114 High mountain E 27 Sta Maria Val Muestair Grigioni 46.600 10.433 1383 Valley bottom 28 Cogolo Pont Trento 46.365 10.689 1190 Valley bottom N-S 29 Passo del Tonale (old) Trento 46.263 10.611 1795 Mid mountain Pass 30 Malga Mare Trento 46.414 10.681 1950 Mid mountain SEFigure 4
Bormio station annual cumulated precipitation: the homogenised series (blue solid line) is compared with the original series (red dashed line). The homogenised period 1930–1959 is clearly visible, when the original series underestimates precipitation.Figure 5
The inner area with stations, numbered following Table4, the smoothed DEM with contour lines, and station MAEs as size of station points.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 p at each DEM cell (x,y) by considering its elevation h:(1)px,y=ax,y+bx,y·hx,y,where a(x,y) and b(x,y) 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 (wiparx,y) 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:(2)wiparx,y=e-Δiparx,y2/cpar,where Δiparx,y is the difference of the parameter values between the station i and the grid cell and cpar 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(x,y) is finally obtained as follows:(3)px,y=a+b·hx,y+kTx,y·ε,where a and b are the global regression coefficients, kTx,y 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 Section2.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.Table 5
Forni supraglacial weather station positions in the last years.
Beginning period End period Latitude N Longitude E October 2005 October 2014 46°23′57.8′′ 10°35′24.7′′ October 2014 November 2015 46°23′52.0′′ 10°35′30.2′′ November 2015 Today 46°23′42.4′′ 10°35′24.2′′Figure 6
Inner area grid points (violet dots) and the locations of AWS1 Forni at Forni Glacier (red dots), superimposed to Google Earth© satellite images.


## 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).Figure 1
Study area domains: grey shaded the inner box area (longitude E 10°10.98′–10°49.14′ and latitude N 46°15.6′–46°37.34′) centred on the Forni Glacier, with the solid black square bordering the outer area (longitude E 8°36′–12°00′ and latitude N 45°00′–47°48′). Black dots indicate the weather stations considered.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 Figure5 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 thehigh mountain environment, whose lower limit roughly corresponds to the upper limit of vegetation [38]; 8 stations are located on slopesat middle altitude(below the treeline) and 6 more stations are located atPasses and plateau. The remaining 14 stations are positioned at thevalley 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 Table1. 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.Table 1
Meteorological networks considered in this work. References for datasets are available at Acknowledgments section.
Network State Region, province ARPA Lombardia, Meteo Italy Lombardia, Sondrio ARPA Lombardia, CMG Italy Lombardia, Sondrio Meteotrentino Italy Trentino Provincia autonoma di Bolzano, Bozen Italy Alto Adige Meteo Svizzera, Meteo Swiss Switzerland GrigioniFigure 2
Availability of weather stations (a) from 1913 to 2015 and station distribution according to the altitude (b), with classes of 100 m. The inner study area is represented in dark green, and the outer area is represented in salmon.
(a) (b)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 (Figure3). 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).Table 2
Reconstruction of absolute precipitation skill for precipitation classes subdivision.
Rain class (mm) Days MAE (mm) RMSE (mm) [0.0–1.0) 326043 0.22 0.96 [1.0–2.5) 32270 1.42 2.10 [2.5–5.0) 27178 2.13 2.86 [5.0–7.5) 17492 2.83 3.69 [7.5–10.0) 10898 3.45 4.54 [10.0–15.0) 14115 3.97 5.26 [15.0–20.0) 8056 4.67 6.22 [20.0–30.0) 7776 5.49 7.40 [30.0–40.0) 3180 6.55 8.87 [40.0–50.0) 1344 7.56 10.31 [50.0–75.0) 984 9.16 12.29 [75.0–100.0) 159 11.92 16.24 [100.0–125.0) 30 13.98 17.48 [125.0–150.0) 7 28.93 35.38 [150.0–200.0) 2 21.67 22.11Figure 3
Scatter plot of simulated versus measured monthly wet days (a) and daily rainfall (b). The skill score is the red line (1 : 1), while blue lines represent the ratios 2 : 1 and 1 : 2.
(a) (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.Table 3
List of homogenised precipitation series and their homogenised periods.
Station name Network Homogenised periods Bormio ARPA Lombardia 1930–1959 Cancano (Valdidentro) ARPA Lombardia 1951–1969 Solda Provincia di Bolzano 1930–1965, 1966–1995 Careser (dam) Meteotrentino 1975–1992 Cogolo Pont (Centrale) Meteotrentino 1933–1940, 1990–1999 Passo del Tonale Meteotrentino 1923–1930, 1959–1971 Mezzana Meteotrentino 1990–2000, 2009–2016 Malga Mare (Centrale) Meteotrentino 1974–2016 Peio Meteotrentino 1930–1950, 1944, 1960–1972Table 4
Inner area station characteristics: the ID column refers to Figure5.
ID Station name Province Lat N Lon E Elevation [m] Positioning Side 1 Grosotto Sondrio 46.286 10.261 613 Valley bottom 2 Premadio centrale Sondrio 46.485 10.351 1305 Valley bottom 3 Trafoi Bolzano 46.550 10.508 1570 Mid mountain E 4 Pezzo Brescia 46.283 10.519 1557 Valley bottom 5 Agumes Prato Bolzano 46.621 10.582 907 Valley bottom 6 Gioveretto dam Bolzano 46.509 10.725 1851 Mid mountain SE 7 Ganda Bolzano 46.550 10.786 1257 Mid mountain NW 8 Lago Verde Bolzano 46.467 10.803 1488 Mid mountain E 9 Rabbi Somrabbi Trento 46.410 10.804 1352 Mid mountain S 10 Santa Caterina Valfurva Sondrio 46.415 10.498 1760 Valley bottom 11 Forni Sondrio 46.421 10.564 2118 Mid mountain Bottom 12 Cancano Sondrio 46.514 10.323 1855 Mid mountain Bottom 13 Bormio Sondrio 46.466 10.371 1225 Valley bottom 14 Solda Bolzano 46.514 10.594 1905 Valley bottom 15 Lasa Bolzano 46.615 10.693 865 Valley bottom 16 Silandro Bolzano 46.622 10.779 698 Valley bottom 17 San Giovanni Martello Bolzano 46.526 10.736 1616 Mid mountain E 18 Pian Palù (dam) Trento 46.337 10.614 1800 Mid mountain NE 19 Careser (dam) Trento 46.423 10.699 2600 High mountain NW 20 Mezzana Trento 46.313 10.796 905 Valley bottom W-E 21 Passo del Tonale Trento 46.263 10.597 1875 Mid mountain Pass 22 Peio Trento 46.364 10.678 1585 Mid mountain SSE 23 Sondalo, Le Prese Sondrio 46.350 10.355 950 Valley bottom 24 Aquilone Sondrio 46.401 10.355 1082 Valley bottom 25 Grosio, Fusino dam Sondrio 46.327 10.246 1220 Mid mountain S 26 Valdisotto Sondrio 46.384 10.330 2114 High mountain E 27 Sta Maria Val Muestair Grigioni 46.600 10.433 1383 Valley bottom 28 Cogolo Pont Trento 46.365 10.689 1190 Valley bottom N-S 29 Passo del Tonale (old) Trento 46.263 10.611 1795 Mid mountain Pass 30 Malga Mare Trento 46.414 10.681 1950 Mid mountain SEFigure 4
Bormio station annual cumulated precipitation: the homogenised series (blue solid line) is compared with the original series (red dashed line). The homogenised period 1930–1959 is clearly visible, when the original series underestimates precipitation.Figure 5
The inner area with stations, numbered following Table4, the smoothed DEM with contour lines, and station MAEs as size of station points.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 p at each DEM cell (x,y) by considering its elevation h:(1)px,y=ax,y+bx,y·hx,y,where a(x,y) and b(x,y) 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 (wiparx,y) 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:(2)wiparx,y=e-Δiparx,y2/cpar,where Δiparx,y is the difference of the parameter values between the station i and the grid cell and cpar 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(x,y) is finally obtained as follows:(3)px,y=a+b·hx,y+kTx,y·ε,where a and b are the global regression coefficients, kTx,y 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 Section2.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.Table 5
Forni supraglacial weather station positions in the last years.
Beginning period End period Latitude N Longitude E October 2005 October 2014 46°23′57.8′′ 10°35′24.7′′ October 2014 November 2015 46°23′52.0′′ 10°35′30.2′′ November 2015 Today 46°23′42.4′′ 10°35′24.2′′Figure 6
Inner area grid points (violet dots) and the locations of AWS1 Forni at Forni Glacier (red dots), superimposed to Google Earth© satellite images.


## 3. Results and Discussion

### 3.1. Climatologies

The monthly climatologies (1961–1990 normals) were evaluated by LWLR for the inner domain reported in Figure1 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.Table 6
LWLR and RK leave-one-out uncertainties. BIAS were evaluated as difference between measured and simulated values.
Month LWLR RK [ mm ] BIAS MAE RMSE BIAS MAE RMSE January −0.48 4.6 5.8 −0.92 7.3 8.8 February −1.52 5.2 6.3 −1.80 7.3 8.5 March −1.20 5.8 7.2 −1.06 8.0 9.5 April −0.94 8.2 10.3 −0.71 11.3 13.8 May −0.81 9.9 12.7 −1.16 9.4 12.4 June −0.78 7.7 9.5 −0.35 7.1 9.8 July −0.16 7.6 9.9 0.13 7.7 9.6 August −1.12 8.3 10.3 −0.57 8.2 10.6 September −0.42 7.7 9.8 −0.58 8.1 11.2 October −0.43 8.3 11.4 −0.91 8.8 11.8 November −1.39 7.2 9.2 −1.12 7.9 10.1 December −0.22 5.1 6.3 −0.46 6.5 7.8Regarding the annual climatological map depicted in Figure7, 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.Figure 7
Annual climatology (1961–1990) evaluated with the LWLR method. Black dots indicate rain gauges; black lines indicate borders between regions. Letters indicate administrative regions: (A) Lombardy, Sondrio province; (B) Switzerland, Grisons canton; (C) Alto Adige Südtirol; (D) Trentino; (E) Lombardy, Brescia province.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 Figure8, 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.Figure 8
LWLR seasons climatologies. Black lines are borders of administrative regions as in Figure6; black dots are weather stations.The same maps of Figures7 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 (Figure11(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 [55–57]. 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.Figure 9
Annual rainfall anomalies. The red solid line represents the short-term trend (11-year Gaussian filter, 3-year standard deviation) and the dashed black line the Theil-Sen trend for the period 1913–2015.Figure 10
Running trends of annual precipitation anomalies. Colours represent the value of trends while pixel size refers to trend significance level. Trend evaluation is performed on windows whose widths range from 20 years up to the end of series starting from each year inx-axis.Figure 11
(a) Scatter plot of monthly simulated LOO precipitation versus the measured precipitations for each station. (b) Scatter plot of stations correlation versus distance.
(a) (b)


### 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 Table4) 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 (Figure11(a)) and we found(4)ysim=0.96xobs+1.36mmR2=0.84.The highest common variance R2=0.96 between the reconstructed and the observed record was found for Rabbi Somrabbi (ID 9 in Table 4); the lowest (R2=0.66) 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:(5)ysim=0.96xobs+1.23mmR2=0.86fall-winterysim=0.96xobs+1.54mmR2=0.80spring-summer.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 Figure12(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.Figure 12
(a) Scatter plot of measured precipitations versus simulated precipitation at grid point 2101. The red line represents the linear regression of data, meanwhile the black line represents the skill score 1 : 1. (b) Scatter plot of monthly model discrepancies versus wind speed. Only discrepancies greater than 20 mm are shown.
(a) (b)A deeper temporal analysis highlighted a high seasonality between simulated and measured precipitation (Figure13); 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.Figure 13
Monthly distribution of discrepancies between simulated and measured precipitation at Forni Glacier for the period 2005–2015. Positive values indicate an underestimation by the model.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 Figure12(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.


## 3.1. Climatologies

The monthly climatologies (1961–1990 normals) were evaluated by LWLR for the inner domain reported in Figure1 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.Table 6
LWLR and RK leave-one-out uncertainties. BIAS were evaluated as difference between measured and simulated values.
Month LWLR RK [ mm ] BIAS MAE RMSE BIAS MAE RMSE January −0.48 4.6 5.8 −0.92 7.3 8.8 February −1.52 5.2 6.3 −1.80 7.3 8.5 March −1.20 5.8 7.2 −1.06 8.0 9.5 April −0.94 8.2 10.3 −0.71 11.3 13.8 May −0.81 9.9 12.7 −1.16 9.4 12.4 June −0.78 7.7 9.5 −0.35 7.1 9.8 July −0.16 7.6 9.9 0.13 7.7 9.6 August −1.12 8.3 10.3 −0.57 8.2 10.6 September −0.42 7.7 9.8 −0.58 8.1 11.2 October −0.43 8.3 11.4 −0.91 8.8 11.8 November −1.39 7.2 9.2 −1.12 7.9 10.1 December −0.22 5.1 6.3 −0.46 6.5 7.8Regarding the annual climatological map depicted in Figure7, 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.Figure 7
Annual climatology (1961–1990) evaluated with the LWLR method. Black dots indicate rain gauges; black lines indicate borders between regions. Letters indicate administrative regions: (A) Lombardy, Sondrio province; (B) Switzerland, Grisons canton; (C) Alto Adige Südtirol; (D) Trentino; (E) Lombardy, Brescia province.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 Figure8, 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.Figure 8
LWLR seasons climatologies. Black lines are borders of administrative regions as in Figure6; black dots are weather stations.The same maps of Figures7 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 (Figure11(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 [55–57]. 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.Figure 9
Annual rainfall anomalies. The red solid line represents the short-term trend (11-year Gaussian filter, 3-year standard deviation) and the dashed black line the Theil-Sen trend for the period 1913–2015.Figure 10
Running trends of annual precipitation anomalies. Colours represent the value of trends while pixel size refers to trend significance level. Trend evaluation is performed on windows whose widths range from 20 years up to the end of series starting from each year inx-axis.Figure 11
(a) Scatter plot of monthly simulated LOO precipitation versus the measured precipitations for each station. (b) Scatter plot of stations correlation versus distance.
(a) (b)


## 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 Table4) 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 (Figure11(a)) and we found(4)ysim=0.96xobs+1.36mmR2=0.84.The highest common variance R2=0.96 between the reconstructed and the observed record was found for Rabbi Somrabbi (ID 9 in Table 4); the lowest (R2=0.66) 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:(5)ysim=0.96xobs+1.23mmR2=0.86fall-winterysim=0.96xobs+1.54mmR2=0.80spring-summer.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 Figure12(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.Figure 12
(a) Scatter plot of measured precipitations versus simulated precipitation at grid point 2101. The red line represents the linear regression of data, meanwhile the black line represents the skill score 1 : 1. (b) Scatter plot of monthly model discrepancies versus wind speed. Only discrepancies greater than 20 mm are shown.
(a) (b)A deeper temporal analysis highlighted a high seasonality between simulated and measured precipitation (Figure13); 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.Figure 13
Monthly distribution of discrepancies between simulated and measured precipitation at Forni Glacier for the period 2005–2015. Positive values indicate an underestimation by the model.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 Figure12(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.



---

*Source: 9123814-2018-03-08.xml*