Application of the MUSLE Model and Potential Effects of Climate Change in a Small Alpine Catchment in Northern Italy

Abstract

The evaluation of sediment yield by water erosion taking into consideration the possible impact of climate change is the object of this work, concerning the use of the Modified Universal Soil Loss Equation (MUSLE) in an Italian case study. This empirical model was implemented in a Geographical Information System, taking into account Alpine hydrology and geomorphological and climate parameters, which are crucial in the analysis of the intensity and variability of sediment yield production processes. The case study is the Guerna Creek basin, a small-sized mountain watershed placed in Lombardy, in the South-Central Alps (Northern Italy). In recent decades it has been hit at the same time by floods and erosive phenomena, showing its hydraulic-hydrological weakness. Three future climate change scenarios from 2041 to 2060, around the middle of this century, were built according to CORDEX data referring to three different Representative Concentration Pathways (RCP 2.6, RCP 4.5, RCP 8.5). The findings showed that in the future climate, the sediment yield at the basin scale might change by 24–44% for a single heavy storm in the middle of the current century.

1. Introduction

Soil erosion processes, which normally contribute to the natural evolution of the landscape, might reach very high intensities during heavy storms and ultimately turn into uncontrolled phenomena leading in the worst cases to desertification [1,2,3]. The importance of the role played by the soil as an ecosystem was already stated in the European Soil Charter [4], then strong emphasis on the soil natural ecosystem services was placed by the publication of the Millennium Ecosystem Assessment [5] and the decision of the United Nations to make 2015 the ‘International Year of Soils’ [6].

Soil erosion by water is the result of water action on an erodible medium (soil) and it mainly depends on the climate, soil, geomorphology, land cover and land use [7,8,9]. At the basin scale, the sediment yield consists of the amount of sediment crossing the basin outlet in a fixed time interval. In several sites worldwide, sediment yield can raise challenging sustainability issues for the management of both the ecosystem and the built infrastructure, and this applies especially to lakes and reservoirs [10,11].

On one side extraordinary meteorological events, on the other side land use change and unsustainable agricultural practices might strongly disturb the natural balance between the soil and the atmosphere. Among the several interactions occurring between the atmosphere and the Earth’s surface, soil erosion by water might be one of those strongly enhanced by global warming, increasing mountain areas’ vulnerability, and requiring the adoption of adaptation and mitigation measures [12,13]. In several sites of the Alpine region, an increase in the frequency of heavy storms might add to an increase in the rainfall peak intensity, leading to a higher amount of sediment yield and urging for more sustainable management practices to be adopted. Moreover, higher values of air temperature might favor, on one side, the soil particle detachment directly, and on the other side, the start and spreading of wildfires, indirectly responsible for higher sediment yields [13,14].

The estimate of sediment yield at the basin scale is though still highly uncertain, because of both the complexity of the physical description of the phenomenon and the lack of observations needed for proper model calibration and validation [15]. As a result, the task is usually accomplished through empirical models, mainly referring to the Universal Soil Loss Equation (USLE) [16], or empirically derived conceptual models, aiming at evaluating the effects of high intensity rainfall events, and the comparison of both different and similar case studies, each adding precious pieces of knowledge to the main scheme [17,18]. Empirical models are usually the simplest and most widely employed to assess soil erosion, particularly at the catchment scale [19]. Physical models, such as WEPP [20,21], have been developed and used mainly to simulate the behavior of farmed fields, but they need very detailed information, both on the soil and on the vegetation, which is hardly available in mountain areas.

An empirical model outlining the event-based effects of surface runoff on sediment yield at the basin scale is the so-called Modified Universal Soil Loss Equation (MUSLE; [22]). As shown by [23], most MUSLE applications were conducted in Asia, especially in Iran, North America and Europe, but no systematic study was done in Italy or the Alpine region.

In this work, MUSLE was applied to estimate the sediment yield during heavy storms in a small Alpine watershed and to investigate how this amount might change in future climate scenarios.

The case study is the Guerna creek watershed, featuring a drained area of about 30 km², located in the Central Southern Alps (Lombardy, Bergamo). Both soil loss by water erosion and sediment transport to the basin outlet, where the Guerna Creek joins the Oglio River outflowing from the Lake Iseo, are challenging local farmer communities and the authority managing Lake Iseo and the Oglio River. The same study area was selected in previous work to evaluate soil erosion by water through the Revised Universal Soil Loss Equation (RUSLE, [24]) and the Erosion Potential Method (EPM, [25]) in the current [26] and in the future climate [14,27]. The novelty of this study lies in the event-based sediment yield estimate. In fact, even over a time period spanning several years, total soil erosion might mainly derive from a few highly effective events. Therefore, soil protection interventions designed on the basis of the estimate of the yearly average soil erosion might not be suitable to cope with the most intense events [28]. Moreover, the climate factor in MUSLE, instead of the one used in RUSLE, aims at the evaluation of the basin sediment yield (rather than of the parcel soil loss), accounting also for solid transport and deposition processes.

2. Materials and Methods

2.1. Study Location: The Guerna Creek Basin

The research was carried out with reference to the Guerna catchment, a small Alpine watershed, 30.9 km² wide, sited in Lombardy, in Bergamo province (North-Central Italy, Figure 1).

The Guerna creek enters the Oglio river from the right and it springs at 1332 m a.s.l., from Monte Foppa. It flows for approximately 12.3 km and, just a few meters downstream the Sarnico dam which controls the level of water in Lake Iseo, it flows into the Oglio river (45.66° N, 9.94° E) [27]. Most of the basin lies in the Municipalities of Viadanica, Adrara S. Martino, and Adrara S. Rocco. Its average slope is around 51% (27°) and the average elevation is approximately 649 m a.s.l. The maximum and minimum elevation are, respectively, 1332 and 185 m a.s.l.; therefore, the range in altitude is 1147 m. Additional characteristics of the Guerna creek and its catchment are listed in

Table 1. Characteristics of the Guerna creek catchment [27,26].

Parameter	Value
Wooded, reforested and forested area	20.5 km2
Urban area	1.83 km2
Grassland	1.64 km2
Bare area 1	0.05 km2
Orchard and vineyard	1.07 km2
Agricultural area	6.00 km2
Longitudinal average slope of the Guerna creek	0.07 m/m
Mean annual precipitation 2	1250.6 mm
Mean annual temperature 2	13.9 °C

¹ Bare area is characterized by a plant cover almost non-existent (for example degraded areas, debris accumulation, rocky outcrops); ² At Sarnico station (years 2008–2011).

. Figure 2 shows, respectively, the land use map and the digital elevation model (DEM) of the study area.

Thematic maps (slope map, lithological map hydrogeological instability map) of the Guerna watershed were produced using a GIS application and they are presented in [26]. In the lithological map limestones, marly limestones, and flint limestones are identified in most of the study area, while gravel, sand, blocks, and ferritization silt are evident at the Guerna mouth. The hydrogeological instability map shows that in much of the watershed, there are no extended erosion effects, even if rill and gully erosion, landslides, debris or terrigenous covering, and alluvial cone can be identified in some spots. These erosion points in the study area were also reported by local authorities and they were noted during field surveys.

2.2. The MUSLE Model and Its Application

The empirical model adopted in this study to estimate sediment yield at the basin scale is the so-called Modified Universal Soil Loss Equation (MUSLE; [22]), which takes into account the runoff effects on soil detachment through the product between the peak discharge and the runoff volume for a single heavy storm.

Ref. [22] used data collected from 18 small US watersheds to evaluate the parameters and the structure of MUSLE. He showed that runoff is highly correlated with sediment yield and therefore the rainfall factor R previously used in the RUSLE can be replaced by a runoff factor Rd [29], accounting for both the climate and the solid transport mitigation effects along the slope. It should be recalled that the RUSLE model is an erosion prediction and conservation planning tool designed to predict the longtime average soil losses; it estimates the mean annual volumes of soil erosion at the basin scale, not the sediment for a single rainfall event. For these reasons, MUSLE and RUSLE models are different and therefore their results cannot be compared. Ref. [22] evaluated the sediment yield (Ys), defined as the amount of sediment crossing the outlet of the basin in a certain interval of time [30], for a single storm using MUSLE:

$${\mathrm{Y}}_{\mathrm{s}}={ \mathrm{R}}_{\mathrm{d}} · \mathrm{K}· \mathrm{LS} · \mathrm{C} · \mathrm{P}$$

(1)

where Y_s (t·ha⁻¹) is the sediment yield, K (t·ha·h·ha⁻¹·MJ⁻¹ mm⁻¹) is the soil erodibility factor, LS (-) is the topographic factor combining slope length L and slope steepness S, C (-) is the cover-management factor, P (-) is the supporting practices factor, and R_d (t·ha⁻¹·(unit of K)⁻¹) is the event runoff factor computed through the following formula:

$${\mathrm{R}}_{\mathrm{d}}=\frac{89.6}{{\mathrm{A}}_{\mathrm{w}}}· {\left(\mathrm{Q}·\mathrm{V}\right)}^{0.56}$$

(2)

where A_w (ha) is the watershed area, Q (m³·s⁻¹) is the peak flow rate of the runoff event, and V (m³) is the volume of the runoff event.

In this research, the sediment yield was estimated for the 4 September 2011 rainfall event at the Sarnico station (see Figure 1); the 2011 event showing the highest recorded peak rainfall intensity, even if it cannot be considered a flash flood event.

The methodology employed in computing K, LS, C and P factors is the same used in [26] to simulate the soil loss by RUSLE in the years 2008–2011 and it is briefly recalled below. The K factor was considered uniform, whereas the LS, P, and C factors were spatially variable in the Guerna catchment. The K factor value is 0.038 t·ha·h·ha⁻¹·MJ⁻¹ mm⁻¹. It was calculated using both equations in [31], which were suggested by [32] and improved in [33], and the outcomes of the particle analysis of soil samples described in [26]. The LS factor was assessed taking into account its spatial variability and adjusting Moore and Burch’s equation [34], by using Nearing’s formula [35]. The LS factor map was built in a GIS environment employing the slope map and the flow accumulation map. The mean, minimum and maximum values of the LS factor are, respectively, 14.09, 0.027 and 41.78. Concerning the C and P factors, the average literature values were used and their spatial variability was considered. The C factor map and the P factor map were created in a GIS environment according to the land use map. The mean, minimum and maximum values of the C factor are, respectively 0.105, 0.002 and 0.45. The mean, minimum and maximum values of the P factor are, respectively 0.98, 0.45 and 1.

For the MUSLE application with reference to the 4 September 2011 rainfall event at the Sarnico station (see Figure 1), the R_d factor was considered uniform over the Guerna catchment. The precipitation amount for each half an hour was provided by Consorzio dell’Oglio [36]. These data were first used to evaluate the surface runoff, using the Soil Conservation Service-Curve Number (SCS-CN) method [37], and then to build the runoff hydrograph at the outlet of the Guerna watershed, using the time–area method [38], from which the runoff volume and the peak flow can be drawn.

2.2.1. Surface Runoff

The first step in building the hydrograph generated by the 4 September 2011 rainfall event, is to estimate the surface runoff rate from rainfall intensity; to achieve this aim, the Soil Conservation Service-Curve Number (SCS-CN) method [37] was employed, which assumes a direct proportionality between soil retention and surface runoff. This method has several positive features and it has already been applied in combination with the MUSLE equation for other case studies, for example in Black Hawk County, USA [39]. Moreover, the SCS-CN model is a well-known and widely applied model in many Mediterranean countries [40,41,42]. The most important benefits of this method are the following: simplicity, predictability, stability, reliance on only one parameter, and responsiveness to major runoff-producing watershed properties (soil type, land use, surface condition, etc.) [43].

In this study, considering that the precipitation fallen in the previous five days was 19 mm and acting in favor of security, the AMC (Antecedent Moisture Condition) was assumed to be of the II class. A table with different CN values for AMC II is given in [37]. In this table, various CN values under AMC II are presented according to the land use and the Hydrological Soil Group (HSG). There are four HSG groups, which are identified as A, B, C, and D on the basis of their minimum infiltration rate and runoff potential. Regarding the Guerna watershed, group D was selected because the results of the particle analysis described in [26] revealed that the soil fraction, as a percentage of the weight passing at 2 mm diameter, is mostly clay and silt.

Table 2. CN values in the Guerna creek basin.

Description	CN (-)
Wooded, reforested and forested area	74
Grassland	71
Agricultural area	83
Orchard and vineyard	83
Urban area	90
Bare area	90

shows the CN values attributed to the different land use areas in the Guerna catchment:

Using values in Table 2, the weighted average value of CN for the watershed is 76.8: this value was used to estimate the surface runoff.

Figure 3 compares the rate of the surface runoff, that was calculated using the SCS-CN method, to the rainfall intensity recorded at the Sarnico station for the 4 September 2011 rainfall event.

Finally, accounting for the high uncertainties in precipitation projections for the future climate, it was deemed useful to also evaluate AMC of the III class to create a more adverse scenario. Usually, most extreme flood events take place after consecutive days of rain, forming AMC III; in this case, the average CN value for the watershed is 88.4. However, it should be noted that the flood event analyzed in this work cannot be considered an extreme event, although it is critical and characterized by a high erosive potential. Moreover, as mentioned above, precipitation fallen in the previous five days was 19 mm. Therefore, AMC II returns a more plausible and still safe scenario. For these reasons, in this work, only AMC II was used to estimate surface runoff and AMC III was considered just to check the worst potential scenario.

2.2.2. Time of Concentration and Time–Area Curve

The second step in the hydrograph building process requires the estimate of the time of concentration in the Guerna watershed and of the time–area curve.

The travel time of surface runoff, or time of concentration, measures the time response of a catchment to a precipitation event. It is the time needed by water to flow from the ‘hydrologically’ farthest point in a basin to its outlet. It is a function of the basin land use, topography and geology [38].

In the technical literature, there are many empirical formulas that find ample application, due to the limited information they need in evaluating the time of concentration on a basin scale in ungauged conditions. Ref. [44] proposed the dominant formula in Italy, which can be used for catchments characterized by drainage areas ranging between 170 and 70,000 km²:

$${\mathrm{t}}_{\mathrm{c}}=\frac{4·\sqrt{\mathrm{A}}+1.5·\mathrm{L}}{0.8·\sqrt{{\mathrm{H}}_{\mathrm{m}}}}$$

(3)

where t_c is the time of concentration (h), A the watershed area (km²), L the length of the main channel (km) and H_m the difference between the mean basin elevation and the outlet elevation (m).

However, the Giandotti formula is not very suited for the Guerna catchment because of its size. The most appropriate formulas for the study area are Tournon [45] and Fattorelli and Marchi [46] ones. The Tournon formula was deduced from observations in north-western Italian (Piedmontese) mountain catchments characterized by drainage areas ranging between 30 and 170 km²:

$${\mathrm{t}}_{\mathrm{c}}=0.396·\frac{\mathrm{L}}{\sqrt{{\mathrm{i}}_{\mathrm{a}}}}· {\left(\frac{\mathrm{A}·\sqrt{{\mathrm{i}}_{\mathrm{a}}}}{{\mathrm{L}}^2·\sqrt{{\mathrm{i}}_{\mathrm{m}}}}\right)}^{0.72}$$

(4)

$$\sqrt{{\mathrm{i}}_{\mathrm{a}}}=\frac{\mathrm{L}}{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}\frac{{\mathrm{L}}_{\mathrm{i}}}{\sqrt{{\mathrm{i}}_{\mathrm{ai}}}}}$$

(5)

where t_c is the time of concentration (h), A the watershed area (km²), L the length of the main channel (km), i_m the basin average slope (m/m) and i_a the main channel average slope (m/m), which is composed of n-stretches of river (each stretch has an average slope i_ai (m/m) and length L_i (m)).

The Fattorelli and Marchi formula was deducted from observations in Alpine basins characterized by drainage areas ranging between 7 and 200 km²:

$${\mathrm{t}}_{\mathrm{c}}=6·\frac{{\mathrm{L}}^{2/3}}{{\mathrm{d}}^{1/3}}$$

(6)

where t_c is the time of concentration (h), L the length of the main channel (km) and d is the elevation range of the main channel (m).

Table 3. Time of concentration for the Guerna watershed by applying different empirical formulas.

Empirical Formula	Time of Concentration (h)
Giandotti	2.21
Tournon	3.15
Fattorelli and Marchi	3.07

shows the time of concentration for the Guerna catchment, resulting from the application of the different empirical formulas mentioned above.

The time of concentration map for the Guerna catchment was built in a GIS environment; the GIS software GRASS was used [47]. Although GRASS GIS works with both raster and vector files, raster data structures were employed because of their analytical power. The elevation raster map resolution adopted to calculate the travel time is 20 × 20 m [48]. The GRASS GIS function used is “r.traveltime” [49,50,51,52], which makes it possible to compute the travel time of surface and channel runoff to the outlet. The travel time from each point inside the basin to the outlet was calculated using the raster GIS functions, which define the flow path and the travel time along the path for each cell. A drainage area, characterized by a threshold value, considers either surface runoff (with overland velocity) or channel runoff (with channel velocity). It was assumed that any upstream area smaller than the threshold value was not able to produce enough runoff to support a channel. Considering the Guerna catchment, 12 cells (0.48 ha) were indicated as channel flow threshold. This choice is in line with the discussion in [26] to estimate flow accumulation value, which was useful to assess the LS factor. Each cell with less upstream area was assumed to be controlled by overland flow.

The runoff velocity is necessary to estimate the travel time through each cell. The GRASS function considers the following elements in determining flow velocities: land use, watershed slope and its location, rainfall excess intensity and channel features. In accordance with [52], the excess rainfall is the total rainfall minus infiltration and other abstractions. The kinematic wave approximation was assumed to derive the runoff velocity for areas characterized by overland flow: the flow depth at equilibrium [53] was used in Manning’s equation to compute the equilibrium overland runoff velocity. An equilibrium discharge [54] for each cell was calculated (Q = Area × Specific discharge) to achieve channel flow velocities and was combined with Manning’s equation. For this study area, the specific discharge used is 5.555 m³/(s∙km²), which corresponds to approximately 19 mm/h intensity of excess rainfall. Considering the precipitation observations recorded at Sarnico station, the maximum rainfall intensity recorded value is 28.8 mm/h, which refers to the study rainfall event (4 September 2011). The travel time through each grid cell was finally determined as the ratio between the flow distance and the flow velocity. The flow distance was defined according to the flow direction, for both surface and channel flow. The flow direction map was built for the Guerna catchment, using the DEM (Figure 2) and the “r.watershed” GRASS function [55]. Considering the watershed structure, the channel width was set to 2 m.

The adopted Manning’s roughness coefficients 0.025 s∙m^1/3 and 0.035 s∙m^1/3, respectively for channels and overland. These values are in line with the evidence provided by photographs and Manning’s values proposed in the literature [56].

In accordance with the GRASS manual, the minimum slope assumed for flat areas is 0.001 m/m. Indeed, the program adopts this value for smaller slopes [49,50].

Figure 4 shows the time of concentration map for the Guerna catchment. The cumulative time of concentration of the entire basin is 3.08 h, in line with Tournon’s and Fattorelli and Marchi’s empirical formulas in Table 3.

The cumulative travel time map (Figure 4) is divided into isochrones which are utilized to produce the time–area curve (Figure 5).

The clear abrupt jump both in Figure 4 (from green to blue) and in Figure 5 (from 2.5 to 3 h) is related to the watershed morphology and to the imposed threshold in the “r.traveltime” function to identify runoff channels. From the analysis of the DEM raster map, it was noted that in correspondence of the abrupt jump there is a significant elevation increase. Indeed, the line of the abrupt jump coincides with the boundary of other sub-basins.

2.2.3. Runoff Hydrograph

The surface runoff event shown in Figure 3 and the time–area curves described in Section 2.2.2 were used to create the runoff hydrograph at the outlet of the Guerna basin. The time–area method [38] was used, based on a convolution of excess rainfall with a time–area curve that represents the progressive area contributions within the watershed in fixed time increments. Pure translation of the surface runoff to the outlet is described using the travel time concept and it results in an outflow hydrograph that ignores storage effects. Using the time–area curve, the translated inflow hydrograph ordinates (Qj) for any selected design hyetograph (that is the rate of surface runoff or excess rainfall in Figure 3) were determined. Each “block” of excess rainfall was assumed to affect the entire Guerna catchment. The simultaneous arrival of the runoff (Q_j) from area A₁, A₂, … for excess rainfall intensity I₁, I₂, … was determined by properly lagging and adding contributions. Therefore, Q_j can be generally expressed as:

$${\mathrm{Q}}_{\mathrm{j}}={\mathrm{I}}_{\mathrm{j}}·{\mathrm{A}}_1+{\mathrm{I}}_{\mathrm{j}-1}·{\mathrm{A}}_2+\cdots +{\mathrm{I}}_1·{\mathrm{A}}_{\mathrm{j}}$$

(7)

where Q_j is the flow hydrograph ordinate (m³/s), I_j is the rainfall excess (or actual runoff) ith-ordinate (m/s), A_j is the time–area curve ordinate (m²) and j is the number of the isochrones contributing to the outlet [38].

Figure 6 shows the simulated hydrograph at the Guerna catchment outlet and the rainfall that caused it.

The peak flow rate (Q = 101.2 m³/s) and the volume of runoff (V = 1086343 m³) of this flood event were used to compute the event runoff factor Rd. The Rd factor obtained by applying Williams’ equation (Equations (2)) is Rd = 922.6 t/(unit K·ha).

2.3. Climate Change Scenarios

In the same way as in [27], different climate change scenarios were analyzed to tentatively examine the effect on sediment yield production by water erosion at the small basin scale. The climate variable used in MUSLE is precipitation.

Coordinated Regional Downscaling Experiment (CORDEX) data were utilized to consider the effect of climate change in the future. CORDEX data contain time series of values of several parameters (precipitation, air temperature, etc.) at specific grid points [57].

The selected Global Climate Models (GCM), which can return projections of how the climate system may change in the future on a scale of about 1000 by 1000 km using mathematical description, was the ICHEC-EC-EARTH [58]. A Regional Climate Models (RCM), which was applied over a limited area and driven by GCM, was necessary to provide information on much higher spatial resolution: RCA4 (the surface processes of the Rossby Centre regional atmospheric climate model) was the chosen RCM [59]. RCM and the dynamic downscaling technique were used to select CORDEX data [57,60]. The domain used in this analysis, which is the region where the regional downscaling takes place, was EUR-11i. The resolution of this domain is 0.125 degrees, both in longitude and in latitude [60].

Therefore, based on the GCM, the RCM and the domain selected, three climate scenarios (or climate projections) were analyzed for the middle of the current century, from 2041 to 2060, according to the Representative Concentration Pathways (RCP) set by IPCC (Intergovernmental Panel on Climate Change) [61]: RCP 2.6, RCP 4.5, and RCP 8.5. In the EURO-CORDEX ensemble, RCP scenarios define pathways of the additional radiative forcing, in W/m², produced by anthropogenic activity (environmental development, socio-economic aspects, and technological and demographic factors) until the end of the current century (the value in 1750 is considered as a reference). RCP 8.5 denotes very high greenhouse gas emission radiative forcing (8.5 W/m²), which continues to increase even after 2100. RCP 4.5 has a radiative forcing value of 4.5 W/m² and it is considered a stabilization scenario, because the force will become stable around the end of this century. RCP 2.6 has a radiative forcing value of 2.6 W/m² and it is a scenario with a significant negative future emission [62].

The observed climate is not synchronized with the CORDEX experiment and therefore the historical simulation had to be examined. In this study, the historical scenario (from 1986 to 2005) in the CORDEX dataset was deemed. Future local climate scenarios were built using observed data at Sarnico (from 2008 to 2011) and comparing historical and future simulations.

2.4. The Impact of Climate Change

As described in [27], by comparing monthly precipitation historical simulations (CORDEX dataset from 1986 to 2005) and future projections (CORDEX dataset from 2041 to 2060 for three different climate scenarios) for the RCM grid point belonging to the Guerna catchment, a correction factor (multiplicative coefficient J) was calculated for each month of the year. The correction factor J for September is the following [27]:

• RCP 2.6 scenario: J coefficient = 0.85
• RCP 4.5 scenario: J coefficient = 0.76
• RCP 8.5 scenario: J coefficient = 1.11

This multiplicative coefficient was then applied, considering each future scenario, to the semi-hourly 4 September 2011 rainfall event in Sarnico station (provided by “Consorzio dell’Oglio” [36]), to project this precipitation into the future and then to evaluate climate change effects. For simplicity and greater clarity, the hourly hyetographs were built in Figure 7, which shows a comparison between different hyetographs, considering or not the effect of climate change:

Hyetographs for RCP 2.6, RCP 4.5 and RCP 8.5 were used to evaluate the hydrographs of surface runoff using the SCS-CN method [37], as described in Section 2.2.1. Considering that the total precipitation amount of the previous five days is 16 mm, 14 mm and 21 mm, respectively for RCP 2.6, RCP 4.5 and RCP 8.5 scenario and acting in favor of security, the AMC (Antecedent Moisture Condition) was again assumed to be in the II class. Here, too, considering climate change impact, it was considered useful to also evaluate AMC of the III class to create a more adverse scenario. Nevertheless, for the reasons already indicated in Section 2.2.1, AMC II describes a more plausible situation and then only this one was applied in this study; AMC III was assessed just for a deeper analysis of the sensitivity of the methodology. Moreover, the increased frequency of short and heavy rainfall events in future scenarios predicted by IPCC [61] suggests that precipitation events might have a lower probability of falling in the III class (AMC III), requesting higher rainfall magnitude on consecutive days previous to the event. CN values attributed to the different land use areas in the Guerna watershed are the same used without considering the impact of climate change (see Table 2); therefore, considering AMC II, the average CN value is 76.8 for all future scenarios and this value was used to create the graphs of surface runoff. Similarly, considering AMC III only for comparison, the average CN value for the watershed is 88.4. The choice not to change CN value in future projections is justified by previous studies [27] in the area, in which it is highlighted that the effects of climate change on land use can be ignored. This conclusion was reached taking into consideration different elements: the comparison between the land use map in the past and the land use map in recent years, the analysis of observed annual precipitation and temperature trends over the years and their comparison considering climate change scenarios. Figure 8 shows a comparison between different graphs of surface runoff, considering or not the effect of climate change:

Then, the graphs of surface runoff were used to create the hydrographs at the outlet of the Guerna watershed for the three future scenarios, using the time–area method [38] described in Section 2.2.2 and the same time–area curve generated without considering the impact of climate change (Figure 5).

Figure 9 illustrates different simulated hydrographs, with and without the impact of climate change:

The peak flow rate and the volume of runoff of this flood event, were used to calculate the event runoff factor Rd for each scenario by implementing Equation (2) (

Table 4. Peak flow rate (Q) and volume (V) and event runoff factor (R_d) for the various scenarios, considering the flood event generated by 4 September 2011 rainfall event and using AMC II to produce surface runoff at Sarnico station.

Scenario	Q (m3/s)	V (m3)	Rd (t/(unit K·ha))
No climate change effect	101.2	1,086,343	922.6
RCP 2.6	76.5	803,094.2	666
RCP 4.5	60.8	628,320.2	510.4
RCP 8.5	120.1	1,308,462	1127

).

Just for comparison, as explained above,

Table 5. Peak flow rate (Q) and volume (V) and event runoff factor (R_d) for the various scenarios, considering the flood event generated by 4 September 2011 rainfall event and using AMC III to produce surface runoff at Sarnico station.

Scenario	Q (m3/s)	V (m3)	Rd (t/(unit K·ha))
No climate change effect	146.6	1,773,341	1494.5
RCP 2.6	119.5	1,413,474	1173.5
RCP 4.5	101.5	1,180,669	968.6
RCP 8.5	166.8	2,044,995	1739.6

shows the same results listed in Table 4, but considering AMC III to calculate surface runoff.

The other factors (K factor, LS factor, C factor, and P factor) in Equation (1) to implement the MUSLE model were not changed compared to values reported in Section 2.1, as they are assumed to be not influenced by climate change. This choice is justified by the fact that the effects of climate change on land use in the Guerna catchment can be neglected, as evidenced in [27].

3. Results

3.1. The MUSLE Modelling of Sediment Yield without Climate Change Effects

On the basis of this analysis and the methods described in Section 2.2, the map of sediment yield (Y_s) as the result of the 4 September 2011 rainfall event (Figure 10) was first built for the Guerna watershed through Equation (1) without considering the climate change. The mean value ( Y_{s_mean}) and the total value (Y_{s_tot}) in the study area are the following:

• Y_{s_mean} = 34.3 t/ha
• Y_{s_tot} = 106,062.4 t

3.2. The Impact of Climate Change According to the MUSLE Model

The sediment yield (Y_s) map was then estimated by solving Equation (1) also considering climate change effects described in Section 2.3 and in accordance with the procedure defined in Section 2.4 (Figure 11). The Guerna basin mean value (Y_{s_mean}) and total value (Y_{s_tot}) of sediment yield, induced by the 4 September 2011 rainfall event, are shown in

Table 6. Sediment yield (MUSLE model) in the Guerna catchment, considering the flood event generated by 4 September 2011 rainfall event and different future scenarios (where AMC II was used to generate surface runoff).

Future Scenario	Ys_mean (t/ha)	Ys_tot (t)
RCP 2.6	24.8	76,593.8
RCP 4.5	19.2	59,401.2
RCP 8.5	42.4	131,171.1

for the three future scenarios.

3.3. Comparison of Scenarios with and without Climate Change and Statistical Treatment

Table 7. Summary of the results: sediment yield generated by the 4 September 2011 rainfall event (MUSLE model) in the Guerna catchment, considering or not the effects of climate change (where AMC II was used to generate surface runoff).

Scenario	Ys_mean (t/ha)	Ys_tot (t)	Ys Variation Compared with “No Climate Change Scenario” (%)
No climate change scenario	34.3	106,062.4	-
RCP 2.6	24.8	76,593.8	−28
RCP 4.5	19.2	59,401.2	−44
RCP 8.5	42.4	131,171.1	+24

summarizes the results achieved for the current and the future climate applying the MUSLE model in the study area and considering the 4 September 2011 rainfall event.

Table 8. Results of the statistical trend of the data obtained by the application of the MUSLE model (4 September 2011 rainfall event) to the Guerna catchment, considering different climate scenarios (where AMC II was used to generate surface runoff).

Scenario	Sediment Yield (Ys) (t/ha) According to MUSLE
	Maximum Value	Mean Value	Standard Deviation
No climate change scenario	635.9	34.3	78.9
RCP 2.6	459	24.8	56.9
RCP 4.5	334.1	19.2	44
RCP 8.5	737.7	42.4	97.1

contains the statistical treatment of the results shown in Figure 10 and Figure 11.

As detailed above, just for comparison,

Table 9. Summary of the results: sediment yield generated by the 4 September 2011 rainfall event (MUSLE model) in the Guerna catchment, considering or not the effects of climate change (where AMC III was used to generate surface runoff).

Scenario	Ys_mean (t/ha)	Ys_tot (t)	Ys Variation Compared with AMC II to Estimate Surface Runoff (%)
No climate change scenario	55.6	171,864.5	+62.0
RCP 2.6	43.6	134,943.6	+76.2
RCP 4.5	36.0	111,381.0	+87.5
RCP 8.5	64.7	200,065.3	+52.5

lists the same results showed in Table 7, but considering AMC III to calculate surface runoff.

4. Discussion

The first step of this work is the implementation of MUSLE [22] in a GIS environment for the Guerna catchment, by solving Equations (1) and (2). The sediment yield was estimated for the 4 September 2011rainfall event at the Sarnico station, and the result obtained is 106062.4 t. To get to this point, rainfall data were used to create the hydrograph at the outlet of the study area, in order to find the volume of runoff (V = 1,086,343 m³) and the peak flow rate (Q = 101.2 m³/s) and to calculate the runoff R_d factor (922.6 t/(unit K·ha)). At this stage, the Soil Conservation Service-Curve Number (SCS-CN) method [37] was used to determine surface runoff from rainfall and, assuming AMC II, 76.8 is the mean value of CN in the study area. AMC III was also considered for comparison purposes and to create a more adverse situation, but corresponding results were not considered because AMC II describes a situation that is more likely to happen in the study area and that already improves security. Moreover, the time–area method was employed to build the runoff hydrograph at the outlet of the Guerna catchment and therefore the time–area curve was generated using the cumulative travel time map built in a GIS environment. The cumulative time of concentration of the entire watershed is 3.08 h and this value is in line with both the Tournon empirical formula and the Fattorelli and Marchi one. The other factors in Equation (1) are the same as reported in [26].

As mentioned above, the scientific literature does not provide useful information to compare the results obtained using MUSLE to other case studies with similar features in the same geographical region. Therefore, examples of MUSLE application in other geographical areas are also considered for comparison purposes, as is the case study in USA proposed by [39], a catchment situated in Black Hawk County, Iowa, with an area of 24.2 km². In this study, the MUSLE model was implemented in a GIS framework and the peak runoff rate was estimated on the basis of both the widely used Soil Conservation Service (SCS) curve number method [37] and the travel time concept. The LS factor was achieved from the DEM, applying the relation of Moore and Burch [34]. The soil erodibility factor K, the cover management factor C and the erosion control practice factor P were derived from Soil Survey Geographic (SSURGO) data and land cover data. The MUSLE equation was finally employed to estimate soil erosion for a typical 24-h rainfall event for a 2-year return period, that is 80.3 mm. Results showed that the total amount of expected sediment yield for the whole basin was 6669 t and the total runoff volume was 765370 m³. In this literature case study, there is much in common with the application of the MUSLE model in the Guerna catchment: the use of Williams’ equation [22] to integrate the model, the usage of a GIS framework, the use of Soil Conservation Service (SCS) curve number method [37], the relation of Moore and Burch [34] to calculate LS factor and the time of concentration to estimate the peak discharge. However, the amount of sediment yield and the total runoff volume are significantly lower than the results obtained in the Guerna catchment without considering the climate change scenario (Table 7). The reason for these findings can be found mainly in the LS factor and the peak discharge, which are both much smaller than those used for the Guerna catchment.

As an example, the application and calibration of the MUSLE model in six small Carpathian watersheds (Poland) [63] are mentioned here. In this case, study, considering the calibration process, the parameters of the runoff factor were estimated using field data. Catchment areas, mainly characterized by forest cover, vary from 32 to 77 km²; depths of rainfall events analyzed range from 4.7 to 53 mm. LS, K, C and P factor values were assessed from the cropping history of the area, using topographic and soil maps and according to literature tables and nomographs. The results show that the product of the peak flow rate Q and the volume V of the runoff event are in the same order of magnitude of value used in the Guerna catchment for the individual storm analyzed. Sediment yield values were estimated in the Polish watershed using both the original form of the Williams’ equation [22] and its calibrated form: in the first case findings are close to the value calculated in the Guerna catchment, but in the second case they are lower.

Other similar examples are the case of Dlubnia and Wisłoka catchments in Poland, analyzed by [64]. The Dublian catchment is mainly covered with forest and the area object of interest is 264 km²; arable land covers most of the Wisloka catchment and the study area is 165 km². In this work, suspended sediment transport was calculated using field data and the results were adopted to evaluate the applicability of the MUSLE model. As in the example already mentioned [63], findings reveal that the product of the peak flow rate Q and the volume V of the runoff event are in the same range of the value used in the Guerna catchment. Sediment yield values were estimated in Dlubnia and Wisłoka catchments using both the original form of the Williams’ equation [22] and its calibrated form: in the former case results are close to the value calculated in the Guerna catchment, but in the latter case they are lower.

Because of the absence of other case studies reported in the literature on the application of MUSLE and the lack of sediment yield observations in the Guerna catchment to validate simulation results without considering climate change scenarios, sediment yield data from a different though similar case studies were considered. The selected similar case study concerns the catchment of the Cordon Creek, a very small mountain basin (5 km²) located in the Dolomites (eastern Italian Alps), where long-term sediment load data were collected. In this case, the mean hillslope gradient is 52% and the annual precipitation is 1100 mm. Climatic conditions in the Cordon Creek watershed are typically alpine: snowpack accumulation and snowmelt runoff prevail from November to May. The mean gradient and the mean width of the main channel are, respectively 13.6% and 5.7 m. On 14 September 1994, an exceptional flash flood event occurred, presenting the maximum water discharge measured (10.4 m³/s) and an hourly averaged bedload intensity much higher than all the other floods. This flood event was characterized by a mean bed load intensity of 225 m³/h [65,66]. Therefore, by considering a specific weight of 2650 kg/m³ [65], the sediment yield value was 596.3 t/h. In the Guerna catchment, the 4 September 2011 study rainfall event produced a peak flow rate of 101.2 m³/s; considering 14.6 h the duration of the hydrograph (Figure 6) and 106062.4 t the total value of sediment yield, the mean hourly sediment yield is 7264.6 t/h. This value is comparable with the observed value in the Cordon Creek watershed, characterized by a drained area and a peak flow rate, respectively 6 and 10 times smaller.

The second step of this work is the application of the MUSLE model in the Guerna catchment as described above, but considering different climate change scenarios (RCP 2.6, RCP 4.5 and RCP 8.5) in the EURO-CORDEX ensemble. A correction factor J was derived for each future scenario, as explained in [27], and it was applied to the 4 September 2011 rainfall event in Sarnico station. Therefore, the precipitation event was projected into the future. The effect of climate change on land use can be ignored, as supported in [27].

Findings revealed that climate change effects on sediment yield by water erosion cannot be overlooked also by the middle of the current century. According to the MUSLE model and RCP 2.6, RCP 4.5 and RCP 8.5, the sediment yield could change by 24–44% on a basin scale. However, the sediment yield is greater than the value assessed without including climate change only for the RCP 8.5 scenario. The reason behind this outcome may be the R_d factor value, which was calculated using the precipitation correction factor J in future scenarios and which is higher than the value assessed without considering climate change only for the RCP 8.5 scenario, as only for this scenario the precipitation is expected to rise in the month of September. However, the analysis reported in [27] shows that the precipitation correction factor J for future scenarios changes every month; therefore, considering another month, precipitation and consequently sediment yield may be higher or lower for each future scenario analyzed. Even though the maximum precipitation correction factors were estimated for February, a rainfall event in September was chosen because it is the month with the highest erosivity index both in the current climate and in the RCP 8.5 future scenario. The specific discharge for the event investigated in the scenario without climate change effects is 3.27 m³/(s·km²), but literature values of specific discharge for common flash flood event are higher, ranging from 9 to 11 m³/(s·km²) [67,68,69]. The aim of this work is though not limited to the analysis of a single flash flood event; rather, it is the study of climate change’s impact on a flood event which causes sediment yield and which reveals critical issues, even if it cannot be considered as an extreme flood event. Especially, the selected rainfall event in the period 2008–2011, which generated the analyzed flood event, is the only one characterized by both the largest total amount of precipitation and the highest erosivity index at the same time within one calendar year (see [27] for more details). Ref. [70] claimed that the combination of variations both in precipitation intensity and precipitation amount may have a stronger impact on soil loss rather than changes in each individual factor. In fact, soil erosion is expected to rise with the increase in precipitation amount and rainfall intensity, on which the rainfall kinetic energy and the erosivity index depend [12,71,72].

The distribution of rainfall in several areas of the same country or continent can either increase or decrease in the future climate [12]; this is consistent with the prediction of IPCC [61] for the global precipitation changes. Other authors [73,74,75] analyzed the long-term impact of climate change on soil erosion by water and sediment yield, evaluating future climate impact to 2100. Ref. [73] examined the future climate change scenarios of soil erosion near Lake Iseo and in Val Camonica, an Alpine valley adjacent to the Guerna basin. The findings of this study reveal that erosion rate in future scenarios can rise or reduce, because of precipitation and temperature which are climatic parameters. Increasing rainfall amounts and intensities usually rise the erosion rate and this is expected to occur also in climate change scenarios, as the MUSLE runoff coefficient depends on the peak flow rate, which is in turn linked to the peak rainfall intensity, and on the runoff volume, which is in turn linked to the precipitation volume. [75] investigated the effects of climate change on runoff and sediment yield in the Aixola catchment (Northern Spain) and they discovered that runoff and sediment yield can increase or decrease. These variations, in line with this work in the Guerna catchment, are strongly correlated with changes in precipitation.

5. Conclusions

The objective of this study was to examine future trends of sediment yield by water erosion in a small Italian Alpine basin, the Guerna watershed, for a single heavy storm, using MUSLE [22] and considering IPCC Representative Concentration Pathways (RCP 2.6, RCP 4.5 and RCP 8.5). Following the MUSLE empirical approach, the runoff factor (R_d factor) was estimated in the current climate and then modified to take into account the impact of climate change in line with the selected future scenarios. The future modifications of the land use were ignored, improving security as evidenced in [27].

The findings revealed that climate change effects on sediment yield production cannot be overlooked: even if in RCP2.6 and RCP4.5 the event sediment yield might decrease up to 28% and 44%, respectively, it might increase up to 24% in RCP 8.5. In fact, sediment yield can increase or decrease in the future, as found in other research papers [73,74,75]. Concerning the validation of the model in the current climate, because of the lack of both other case studies reported in the literature and sediment yield observations in the Guerna basin, a similar case study was considered for comparison purposes [65,66], supporting the reliability of the results found in this work.

This paper presents a first challenge to estimate sediment yield in the Guerna basin, also considering the potential impact of climate change and using an empirical model. It may be used as a departure point to better understand hydrological impacts in the study area and provides precious insights to local stakeholders and policymakers. However, field measurements for the model’s validation and the expected land use changes to build future scenarios could improve the accuracy and consistency of the whole methodology and its implementation. In spite of the current limits, this work contributes to the current literature on the effects of soil erosion by water, also considering the effects of climate change, which is a very crucial issue today.