Modelling seasonal variation of gully erosion at the catchment scale

Geo‐hydrological phenomena, including gullies, contribute significantly to soil erosion and land degradation. To address this, proper management of basins and hillslopes should consider the mechanism, timing, and location of gully development and how gullies interact with other hillslope processes. Yet, conventional modelling techniques for such processes are rare, frequently being limited to applications of only single processes and typically requiring high‐resolution input data. Further, existing tools for characterizing basins and hillslopes tend to be based on static descriptions of geo‐environmental conditions, and thus are not effective for modelling changes such as the seasonal triggering conditions of gully phenomena over time. This study proposes a method to integrate open remote sensing data (Sentinel‐2) and an existing modelling tool (LANDPLANER) using simplified input data to better predict and forecast gullies’ spatial and temporal occurrence. The study investigates the seasonal conditions responsible for the triggering of gullies at the catchment scale using different erosion modelling schema in the Toscana region of Central Italy. Geomorphological gully inventory data were collected and used as benchmarks to test the proposed approach. The results show that the occurrence of gully erosion in the studied region changes seasonally, and the proposed method was able to effectively discriminate spatial and temporal differences of the gully phenomena. The proposed method can be applied to similar regions worldwide, allowing for the investigation of gully erosion over time, even in places with limited data availability.


| INTRODUCTION
Hillslope geo-hydrological processes triggered by rainfall are widespread and cause severe damage to structures, infrastructures, and population. In the Mediterranean region, such processes are widely distributed and severely impact the stability of slopes (Taguas et al., 2015).
In Italy, intense rainfall events triggering landslide phenomena and erosion processes occur in many regions, causing damage and human losses (Salvati et al., 2010(Salvati et al., , 2014. Water erosion is a relevant threat for a large part of Italy; approximately 77% of the territory is affected by this natural process (Gazzolo & Bassi, 1961). Water erosion mainly occurs on steep slopes in places subjected to intense mechanical cultivation, where conservation measures such as terraces and waterways have been removed (Chisci, 1986;Hauge, 1977), as well as along hillslopes where erodible subsoil materials are exposed (Torri et al., 2006). Even the notable high rate of landslide mobilization is the result of agricultural and land use practices that favour erosion (Torri et al., 2006). In addition, soil erosion is the foremost cause of soil degradation worldwide (Lal, 2001) and is capable of severely affecting water quality and ecological health and negatively altering terrestrial and aquatic habitats (Bennett & Wells, 2018). Gully erosion is a geo-hydrological process that strongly impacts slope dynamics, thereby contributing significantly to soil erosion and influencing other hydrologic phenomena, as well as worsening hillslope environmental conditions (Bennett & Wells, 2018). Specifically, a gully is a hillslope incision in which runoff water accumulates, removing soil material even at considerable depths (Poesen et al., 2003). Moreover, gully channels are typically deep enough (>0.5 m) to interfere with normal tillage operations (Bennett & Wells, 2018). According to the literature, a gully is generally typefied as having a cross-sectional area greater than one square foot, approximately equal to 929 cm 2 (Hauge, 1977).
Gullies are commonly generated by surface runoff and tend to be a persistent feature of the landscape that cannot be remediated by tillage (Gilley, 2005). Moreover, gullies contribute significantly to the total soil erosion amount within a catchment and affect the catchment's morphological shaping (Rossi, 2014;Torri, Poesen, et al., 2018).
These relevant scientific and social issues highlight the importance of providing more effective modelling approaches to predict the spatial and temporal occurrence of gullies (Rossi, 2014).
Changes in land use, such as deforestation and the substitution of forests by cropland or meadows, may affect the occurrence of erosion phenomena (Torri, Poesen, et al., 2018). Generally, erosion phenomena are triggered by long-term, complex anthopogenic processes acting over centuries, such as inappropriate cultivation or irrigation practices, over-grazing, and over-cropping, but also by recent practices such as road building and urbanization (Poesen et al., 2003;Valentin et al., 2005). In one study, Poesen et al. (2003) determined that fluctuations of total soil loss caused by gully erosion can be related to changing land use conditions during a given period, with the amount of soil loss by gullying depending on the time span considered. Gully erosion is also known to be controlled by the magnitude and frequency of rain events (i.e. directly influencing the amount of runoff on hillslopes), which characterize a given climate and weather regime. Vanmaercke et al. (2016) analysed the scientific literature and found that, for a given region, the rainy day normal (RDN, a climate index accounting for the long-term average annual rainfall depth divided by the average number of rainy days) is significantly correlated with the gully head retreat rate. As such, any change in the rainfall regime (due to climate changes) for a specific area will likely affect gully erosion. For example, continental areas may be considerably dominated by gully erosion from concentrated snowmelt runoff in the spring, then by sheet and rill erosion from thunderstorms in the summer (Poesen et al., 2003). Worldwide, similar examples illustrate that land use change is expected to have a greater impact on gully erosion than climate change (Valentin et al., 2005), yet this is in conflict with the findings of Vanmaercke et al. (2016), who were unable to identify significant correlations between gully retreat and land use or soil type. Overall, the factors of climate, land use, and land cover cannot be considered independently, yet at present their combined interactions and joint contribution to gully occurrence are still poorly understood. Consequently, from this perspective, it is important to consider seasonal changes that can intensify the susceptibility of a land type to erosion (e.g. the co-occurrence of seasonal maximum precipitation and little soil vegetation cover).
In this study, we investigate different seasonal conditions which have been identified as relevant for the triggering of gully erosion at the catchment scale (Bogen et al., 1994;Casalí et al., 1999;Monsieurs et al., 2015;Tufekcioglu, 2018). Related studies have highlighted how land use, land cover, and meteorologic and climatic seasonal conditions may affect the triggering and development of gullies, their shape, and the related channel geometry. We therefore propose a method to integrate remote sensing data from the Sentinel-2 satellite with an existing modelling tool, LANDPLANER (Landscape, Plants, Landslide, and Erosion), which describes the dynamic response of slopes under different changing scenarios, to better predict and forecast gully spatial and temporal occurrence. This method aims to solve a critical issue in erosion modelling in which models commonly assume an inadequate averaged (i.e. static) temporal description of the geo-hydrological conditions controlling rainfall/runoff repartition at the surface, and in turn erosion. The proposed method considers more dynamic (i.e. seasonal) descriptions of such conditions based on remote-sensed data with a high revisiting time, and uses them to feed the LANDPLANER model (Rossi, 2014), thereby enabling the model to dynamically adjust the hydrological runoff response of hillslopes and better evaluate their capability to trigger erosion. Moreover, the modelling tool was selected for its advantage of requiring only limited input data.
The proposed modelling framework is applied to a portion of the Freddana Torrent catchment in the Toscana region of Central Italy, in which rainfall-induced geo-hydrological processes (i.e. gullies but also landslides) occur frequently and substantially impact structures and infrastructures. To achieve this, we have: (i) configured LANDPLANER with dual-approach erosion modelling schema; (ii) used different satellite images (eight Sentinel-2 images acquired within a period of 2 years at two images per season) for dynamically characterizing the study area conditions; and (iii) collected information on past gully erosion occurrence (i.e. geomorphological gully inventory) as a benchmark to test the proposed approach.
The rest of the paper is outlined as follows. The next section describes the study area and the available gully inventory data. The third section outlines the entire modelling framework, including methodologies, used to characterize and improve the spatial and temporal variation of the model input data used in the analyses.
The fourth section presents the results. Finally, discussion and conclusions are presented in the last two sections, respectively, and indicate the advantages and limitations of the proposed framework to model seasonal variation of gully erosion and related phenomena at the catchment scale.

| STUDY AREA
The study area (7.1 km 2 , Figure 1 According to the Kӧppen-Geiger classification (Peel et al., 2007), the climate is classified as Csa (hot-summer Mediterranean), which is characterized by hot and dry summers with at least one winter month with more than three times the rainfall amount in the driest summer month, and with the average temperature of the hottest summer month greater than 22 C. The climatic characteristics and rainfall regime make this area among the most vulnerable to geo-hydrological processes in the Toscana region. Indeed, geo-hydrological processes occur frequently in this catchment, mainly during intense rainfall events, including two events recorded in July 2014 and December 2017, respectively (Autorità di Bacino Pilota del Fiume Serchio, 2014a). Erosion, and particularly gully erosion, is a relatively common process, that together with landslide activity impacts the hillslopes, structures, and infrastructures therein.
The closest thermometric and rain gauge stations are located, respectively, at $1 km northeast (Aquilea thermometer, E 1619260 N 4859827, 134 m a.s.l.) and 500 m south of the study area (Mutigliano rain gauge, E 1620624 N 4863398, 33 m a.s.l.) (Figure 1a) The geology belongs mostly to the Tuscan Domain (Tuscan Nappe Unit), while formations of the Ligurian Domain also outcrop in the area, and quaternary slope and alluvial deposits characterize the transition between the hills and Freddana Torrent plain (Nardi et al., 1987). In order of extent, the main geological formations in the  calcilutite and radiolarite sequence with argillites (Cerrina Feroni et al., 2017;Puccinelli et al., 2015Puccinelli et al., , 2016 showing that the soil fertility is low (Sanchez et al., 1982) for most of the study area, thereby sustaining a limited range of crops and conservation practices (Table 1). In the vicinity of the torrent banks, the soils allow for a wider range of crops and/or conservation practices. The Hydrologic Soil Group (HSG) classifications (USDA and NRCS, 2015 -Chapter 7) are mostly B (i.e. soils with 10-20% of clay and 50-90% of sand with a loamy sand or sandy loam texture), which is characterized by soils with relatively low surface runoff potential, and some C, which is characterized by soils with relatively high surface runoff potential.
The land use and land cover (LULC) map of the Toscana region is based on the CORINE Land Cover (CLC) Level III classification schema (Bossard et al., 2000;European Environment Agency, 1995) ( Figures 1e and f). Note that the class called '210: irrigated and nonirrigated crops' mentioned in the subsequent discussion is a newly In the analyses, we initially considered both the 2007 and 2013 LULC maps, but the differences among the two were few. Therefore, it was decided to only show results obtained using the more recent 2013 LULC map. Figure 1 and Table 1 show that the study area is mainly covered by forest (classes 311, 312, and 313), followed by permanent crop   (Beer & Johnson, 1965), dynamic, static, deterministic, and stochastic approaches (Sidorchuk, 1998(Sidorchuk, , 1999(Sidorchuk, , 2005(Sidorchuk, , 2006. Some are based on semiempirical physical models (De Ploey, 1992;Flanagan et al., 2001;Foster, 1986;Tucker et al., 2001;Willgoose, 2005). The more complex models use alternate parameters that are difficult to quantify or locate details of in the literature. Using the gully head topographic threshold (in this study denoted as esp) is common for the empirical approaches (Montgomery & Dietrich, 1994;Poesen et al., 2011;Vandekerckhove et al., 2000). Such an approach has recently been refined by Torri and Poesen (2014) and Torri, Poesen, et al. (2018).
Their models aim to predict the locations at which the gully processes cannot continue expanding upslope, that is, identifying the head of the gully, under specific hillslope conditions and a given rainstorm intensity. Such models are focused on the definition of the critical topographic conditions, expressed by the local gradient(s) of slope and drainage area (A), that control the development and position of the gully heads in various landscapes. Torri and Poesen (2014) suggested that the threshold coefficient k (the resistance of the site to the development of the gully head) depends on land use. They also proposed a new relationship (Torri & Poesen, 2014;Torri, Poesen, et al., 2018) for describing overland flow-induced gully heads, extracted from datasets collected in various parts of the world. A key factor of that proposed relationship is that for a given value of the threshold exponent b, the coefficient k can be predicted using the characteristics of soil and veg- predict landslides and erosion processes by simpler means and requiring only basic input data, such as at least one digital elevation model (DEM) and one land cover map to derive runoff parameters based on the CN method (Rossi, 2014;USDA and NRCS, 2015). The model can be used to simulate water repartition at the soil surface (i.e. the runoff and infiltration repartition). The spatial resolution of the model analysis should match the DEM resolution, in this study corresponding to 10 m. Moreover, the model is well suited to describe the dynamic response of hillslopes (i.e. hydrological and geo-hydrological processes variation) under possible topographic, climatic, and LULC changes.
A key input for the model is the CN map for the catchment area, which is generally defined based on a land use and a soil type map (USDA and NRCS, 2015). To better account for the spatial and temporal variation of land use, we have also proposed a methodology to derive CN maps starting from (i) the land use categories classified according to the CORINE Land Cover scheme (Bossard et al., 2000) and (ii) the variation of normalized difference vegetation index (NDVI) values within each category. To achieve this, we used European Space Agency (ESA) Sentinel-2 images for deriving NDVI values, mainly due to their resolution and high revisiting time (Drusch et al., 2012), which allows for a more accurate spatial and temporal characterization of land use changes in a given area at the seasonal scale. It is important to note that the ability to correctly characterize the variation of these model input data (meteorological and CN) in space and time is fundamental to the prediction of gully formation.
In the following sections, we introduce the proposed model and its input data and describe our approach for improving the spatial and temporal characterization of the input CN map, with the main objective being to strengthen the ability of the model to predict runoff and erosion phenomena (gully erosion).
The flowchart in Figure

| Model description
The open-source modelling tool LANDPLANER is mainly designed to describe the dynamic response of slopes under different changing scenarios, including meteorological factors, vegetation, or land use and slope morphology. The model is raster-based and is able to estimate the effects of given rainfall characteristics on the triggering of landslides and erosion processes (Rossi, 2014). The model encompasses three components: (i) a hydrological model component for the water repartition that explicitly considers the rainfall (P), runoff (Q), infiltration (I), exfiltration (Ex), evaporation (Ev), and transpiration (Tr); (ii) a twofold erosion modelling schema to predict locations where erosion processes may occur; and (iii) a landslide model component to predict where landslide phenomena may occur. The erosion model estimations can either be static, based on a morphologic threshold approach which is mainly used for gully head prediction (Torri & Poesen, 2014), or dynamic, based on the use of an erosion index (in this study denoted as e) (Rossi, 2014) built on the stream power concept, which considers direct runoff (calculated by a specific rainfall/runoff routing procedure), local slope, and soil conditions expressed by the abstraction coefficient S 0.05 that is derived directly from the CN value (equation 3-3 in Rossi, 2014). The site index, S, varies from 0 to ∞ and is a measure of the catchment hydrological response (USDA and NRCS, 2015). The erosion index, as opposed to the topographic threshold, does not solely consider gully erosion, but since the erosion index is directly calculated from runoff, it may potentially account for the different types of erosion processes, including sheet and rill erosion.
In this study, the model utilized simplified input data, namely one DEM (with 10 Â 10 m resolution) to derive different morphometric variables (i.e. accumulation, drainage direction, and slope) ( Figure 3), one LULC map, and one soil map to derive first the CN, and then the S 0.05 maps, which were then used by the model to simulate water repartition at the soil surface.
The rainfall-runoff model is based on the CN method (USDA and NRCS, 2015), thereby allowing for the estimation of runoff based on land use and soil condition data. The approach is able to quantify the runoff produced by a specific rainfall event and can be effectively applied to estimate runoff from daily or sub-daily rainfall data. For this work, the most useful characteristic of LANDPLANER is its capability to use both meteorological input data (i.e. rainfall and temperature) and CN values, as these variables directly express the soil and land use conditions. These variables together allow for reproducing the temporal and seasonal variation of the conditions in the study area which are triggering gully erosion phenomena.

| Dynamical input data derivation
The CN map is one element of the key input data, which strongly con- Following this approach, we were able to identify the spatial variation of CN for each LULC/HSG combination in the study area and used that to define the three CN scenarios reported in the maps in Considering that in the study area, the sea wind either does not blow or only blows with a relatively low intensity (Vittorini, 1972), the type of aerosol was determined as 'rural'. The 'Mid_Latitude' parameter is related to the seasonal atmospheric conditions, hence the selection is based on the date and season of the remote sensing images collection T A B L E 2 Correspondence between CLC classes in the study area defined in was applied for remote sensing images from October to March and 'SUMMER' for those from April to September (Figure 2). In addition, to restrict the possible effects of the planetary boundary layer, that is, the portion of the troposphere influenced by the interaction with the earth surface (i.e. topography) (Stull, 1988), we enabled Sen2Cor is strongly affected by changes of surface optical properties (Baret & Guyot, 1991), can be computed as follows (Nemani & Running, 1989;Rouse, 1974): where R NIR and R RED are reflectances in the near-infrared (0.785-0.900 μm) and red visible (0.698-0.785 μm) wavelengths, respectively (Gamon et al., 1992), and B is the corresponding Sentinel-2 MSI band (Clerici et al., 2017). Standard NDVI values range from À1 to +1 (Myneni et al., 1995;Pettorelli et al., 2005), with three correspondence levels as follows: (i)  To verify the NDVI calculation procedure, we analysed the NDVI value distribution in urban areas (corresponding to the CLC class 112) within and adjacent to the study area, using violin plots ( Figure 5).
Good values were obtained for the study area, that is, values equal or close to 0 (i.e. NDVI values were around 0 and median values were close to 0), as shown in Figure 5. The residual variation of NDVI around 0 can be attributed to the limited LULC pixel variability within the urban-area polygons.
It is known that NDVI temporal profiles efficiently capture the photo-synthetically active vegetation behaviour (Benedetti & Rossini, 1993), and thus can be directly linked with plant phenology variations (Benedetti et al., 1991). This relationship can be observed by comparing Figure 1f with the NDVI maps in Figure 6.
Based on this, it can be stated that NDVI values are generally inversely correlated with CN, which takes its maximum value for impervious soils or bare surfaces. This is a key assumption of the

| Experimental setup and analysis
The   Table 3 we considered Sentinel-2 images relative to different seasons; one for spring, three for summer, two for autumn, and two for winter.

| Model evaluation
In the literature, different metrics, such as receiver operating characteristic (ROC) curves (Kariminejad et al., 2020) or other statistical aggregated metrics, have been used for gully model evaluations. Nevertheless, such metrics do not allow for the estimation of the eroded area length as correctly predicted by the model, which is an important factor for linear erosion features, such as the gullies.
Here, to evaluate and verify the proposed methodology, we compared the modelling output with the gully erosion phenomena mapped as linear features in the study area (Figure 1c). We used two metrics, the first investigating only the length of the gullies correctly predicted by the model and the second investigating the overall model performance on the entire study area. For all simulated scenarios in Table 3 for each gully mapped, we evaluated: (i) the total pixels of the gully (i.e. denoted as the variable 'Gully cells'); (ii) whether the gully was entirely predicted ('Yes') or not ('No') (i.e. 'Entire'); (iii) the number of pixels of the gully that were correctly predicted (i.e. 'Cells predicted'); and (iv) the percentage of the gully that was correctly predicted (i.e. '%', calculated as: sum('Cells predicted')/sum('Gully cells') Â 100).
Finally, we counted the number of gullies with a percentage of predicted cells greater than 0% (i.e. meaning at least one cell was predicted correctly), greater than 50% (i.e. at least half of the gully cells were predicted correctly), and equal to 100% (i.e. all gully cells were predicted correctly). These values were then used to derive different binary classifiers, thereby quantifying the model spatial classification performances (Fawcett, 2006) on the entire study area considering the prone or not  (Kariminejad et al., 2020;Meyer et al., 2008;Rossi et al., 2010) were also produced to graphically illustrate the contingency table analysis outputs. While ACC and TR_AVG are expected to quantify the model spatial classification performances regardless of the distinction between prone or not prone gully areas, the TPR (i.e. the proportion of pixels with observed gullies correctly predicted) and TNR (i.e. the proportion of pixels without observed gullies correctly predicted) quantify them separately. Therefore, the higher these binary classifiers with respect to 0.5, the better the performance of the model, with 1 being a perfect prediction. In addition, a good classification model should jointly maximize TPR and TNR to guarantee the same capability of predicting both the prone and not prone gully areas (Kariminejad et al., 2020). These analyses were performed for each scenario, calculating the e for different rainfall scenarios.
Similarly, the contingency tables and the associated binary classifiers were also derived for the scenarios estimating erosion through the esp modelling approach. For this, the binary maps restituted by the model were directly compared with the rasterized binary gully inventory map.

| RESULTS
For all simulated scenarios, LANDPLANER restituted the erosion estimations, specifically the e maps for different rainfall scenarios (10, 40, 70, and 100 mm) and maps of the possible location of gullies based on the esp maps. The cumulative computation time for all analyses of the study area was approximately 12 h using a basic PC desktop configuration. It should be noted that the extension of the analyses to larger study areas is possible, even with a sensible reduction of the computational times, if using a more powerful hardware configuration or dedicated calculus infrastructure.
F I G U R E 9 Results of LANDPLANER simulations in terms of erosion index (e) for a daily rainfall of 100 mm. a-c represent CN values derived only by LULC maps (Table 2) and d-l represent CN values derived using Sentinel-2 images with (d) Equation 3 and (e-l) Equation 4. See Table 3 for scenario references [Color figure can be viewed at wileyonlinelibrary.com]

| Seasonal analysis
The results of the seasonal occurrence of gully erosion using Sentinel-2 images (Table 3)  For each map, we calculated the evaluation metrics described previously, using the gully mapped in Figure 1c as a benchmark. Table 4 summarizes the values of the metrics evaluating the length of the observed gullies that were correctly predicted.

| Spatial classification performance
The binary classifier metrics, applied to the e model outputs and obtained for different e t values, are summarized in Table S1 in the Supplementary Material. Specifically, Table 5 reports an extract of this table for e t = 0.1 and event rainfall equal to 100 mm, showing the best overall classification perfomances. The fourfold plots in Figure 11 provide graphical illustration of the values of the contingency table reported in Table 5 and help interpret the evaluation output. In Figure 11, the quarter circle areas in the upper parts of the plots are proportional to the TN (blue) and FP (red) frequencies reported in the contingency table; the larger the blue area with respect to its horizontal red counterpart, the larger the TNR value. Similarly, the quarter circle areas in the lower part of the plots are proportional to the TP (blue) and FN (red) frequencies, and hence the larger the blue area is with respect to its horizontal red counterpart, the larger the TPR. When both the blue areas are larger than their horizontal red counterparts, the model shows good classification performance (TNR > 0.5 and TPR > 0.5).
The spatial classification performance metrics were also calculated for the esp model outputs and the relative results are summarized in Table 6.
F I G U R E 1 0 Results of LANDPLANER simulations in terms of topographic threshold (esp). a-c represent CN values derived only by LULC maps (Table 2) and d-l represent CN values derived using Sentinel-2 images with (d) Equation 3 and (e-l) Equation 4. See Table 3 for scenario references [Color figure can be viewed at wileyonlinelibrary. com]

| Percentage variation maps
To better highlight the differences between the e results, obtained with CN values derived using Sentinel-2 images (i.e. changing seasonally) and static CN maps derived only from LULC maps (i.e. not changing seasonally), we produced percentage variation maps. In this process, once a benchmark map is assumed (i.e. RCNMEAN here was considered as the benchmark, x benchmark ), the percentage variation of a scenario map (x scenario ) is calculated using Equation 5. Figure 12 shows the calculated values for a selected set of scenarios.
The comparison of the e estimations in terms of percentage variation with respect to the benchmark scenario RCNMEAN effectively reveals the increase or decrease of gully erosion in the area. The maps show that the e obtained from the NDVI-based scenarios in the summer images has lower values than those in the RCNMEAN scenario. In contrast, the e values for the winter images match better with the spatial gully occurrence.

| DISCUSSION
This study has proposed a method to better account for the spatial and temporal distribution of triggering geo-environmental conditions of gully occurrence through the integration of remote sensing Sentinel-2 data and LANDPLANER. To achieve this, we modelled gullies using a twofold erosion modelling schema under different static and dynamic CN input data scenarios and provided different criteria and metrics to evaluate the results.
The first erosion modelling approach, based on the gully head esp concept, did not require calibration and its spatial prediction variability depended primarily on morphometric parameters (e.g. local slope T A B L E 4 Evaluation metrics comparing observed and modelled data for the different scenarios. Scenarios based on CN derived from NDVI data are based on Sentinel-2 remote sensing data. '0' stands for the number of gullies with a percentage of predicted cells greater than 0 (i.e. at least one cell was predicted correctly); '50' stands for the number of gullies with a percentage of predicted cells greater than 50% (i.e. at least half of the gully cells were predicted correctly); 'ENT' stands for the number of gullies with a percentage of predicted cells equal to 100 (i.e. all the gully cells were predicted correctly). In addition, '% TOT' stands for the percentage of the gully correctly predicted (i.e. variable '%' calculated as: sum('Cells predicted')/sum('Gully cells') Â 100)). Bars express graphically the metric differences T A B L E 5 Extract of the contingency table values and related binary classifiers obtained for the erosion index output maps for an event rainfall amount of 100 mm and an e t threshold of 0.1. Full evaluation results are shown in Table S1 in the Supplementary Material. Acronyms as before. The geomorphological gully inventory map was used as a benchmark. * Denotes the scenarios with better spatial classification performance (i.e. jointly maximizing TNR and TPR above 0.5) which was mainly to compare the spatial and temporal erosion patterns for different scenarios, a specific calibration could be avoided and the calibration parameters desumed from previous modelling applications (Rossi, 2014), yet could also quantify the proposed evaluation metrics for different e t values. The evaulation results showed that a value of e t = 0.1 was appropriate for the analysis performed in this study, in that this threshold value provided the best modelling performance ( Table S1 in (Table 4 and  Similarly, the total number of gullies partially and/or entirely predicted by the model using e, increased with the amount of rainfall. This is because rainfall has the primary effect of increasing the amount of runoff (equation 3-4 in Rossi, 2014), which is directly related to the final e value (equation 3-38 in Rossi, 2014). This occurred at the expense of an increased number of false positives (Table S1 in the Supplementary Material, Table 5 and Figure 11), since the models for larger CN and rainfall values tended to predict a larger area affected by gullies (Figures 9 and 10 and Figures S1-S3), with variation of the e values being greater than 100% ( Figure 12).
The results of the RCNMEAN scenario show a higher number of predicted gullies with respect to RCNNDVIMEAN derived using Sentinel-2 remote sensing images ( The analysis of the spatial model classification performances ( Figure 11 and winter (RCNNDVI20161218, RCNNDVI20171220) show on average a higher number of correctly predicted gullies (Table 4) along with improved spatial classification performances ( Figure 11 and Table 5) with respect to RCNMEAN, which should be considered as an average yearly benchmark. In the study area, during these two seasons and  Figure 6) and the higher CN values (Figure 8).  (Table 4) and lower spatial classification performance (Table 5 and Figure 11). Such seasonal differences, which can be observed directly in the NDVI and CN input data (Figures 6   and 8) and in the modelling outputs (Figures 9, 10, and 12), reflect the conditions less prone to the gullies occurrence during spring and summer periods, which in the study area are characterized by the maximum vegetative development for all the LULC classes. In addition, for these two seasons, the maximum daily rainfall is lower with respect to autumn and winter, even if relatively short, intense rainfall events occur.
The results obtained using the esp approach were commonly between those obtained using e for event rainfall values 70 and 100 mm, with the greater differences obtained mainly during autumn and winter. This justifies the use of the esp approach for a rapid yet effective evaluation of gully occurrence in a given area. It should be noted that the use of esp, commonly assuming a static CN input, within the proposed framework maintains the capability to describe seasonal variation of the gully occurrence.
We acknowledge that the evaluations performed in this study, although in favour of the proposed seasonal dynamics-based methodology, are in some cases only slightly better than the corresponding averaged counterparts. We attribute this to the dependence of the results on the quality and completeness of the geomorphological gully inventory used as a benchmark, which certainly can be considered as an underestimation of the true values. In fact, neither the visual interpretation of stereoscopic aerial images nor the analysis of the hillshade at 1 m resolution derived from two LiDAR surveys, even if additionally aided by field surveys, could fully capture all gully occurrences. Such underestimation may also be due to the gully activity normally realized in the study area within the first months after gully occurrence, specifically the ephemeral gullies and/or smaller gullies infilling and levelling. In this respect, future additional efforts should attempt to build multi-temporal and event-based gully inventory maps, through continuous remote sensing and field-based study area monitoring.
Overall, quantifying the impact of the seasonal variability on erosion for each LULC class is no trivial task; however, the proposed methodology has proven effective to account for these effects, as further supported by the summary statistics shown in Figure S17 in the Supplementary Material. We believe that the proposed modelling approach, combined with timely field observations and characterizations, should help to improve our understanding of seasonal gully occurrence.
The results of this study encourage the use of CN for characterizing such diversified seasonal conditions and show the effectiveness of CN in controlling the repartition of rainfall among the different water balance components, which in turn affect erosion. Nevertheless, in the proposed approach for deriving dynamical input CN data, it is relevant to select satellite images with appropriate timing for better representation of the study area conditions controlling gully occurrence.
For this purpose, using freely available Sentinel-2 data, which have a high revisiting time (Drusch et al., 2012) allows for the possibility to select the most appropriate images based on seasonal timing. Potentially, Sentinel-2 data may also allow for the possibility to set up simulations considering up to decadal LULC variation conditioning of the geo-hydrological processes and more specifically, gully occurrence.
Additionally, the spatial resolution of Sentinel-2 data, which is approximately equal to 10 m, together with the proposed methodology of this study to convert NDVI to CN, allow for the possibility to better represent the LULC spatial variability, in that normally the CN from literature is done at lower spatial resolutions, as it is based on the CN variation within LULC patches (polygons) in place of single pixels.
As highlighted by Vanmaercke et al. (2021), a variety of processbased and empirical models have been proposed to predict gully occurrence, expansion, and contribution to sediment yield, and also to estimate gully impacts. They suggest that there is the need for process-oriented model approaches to investigate the factors and mechanisms driving gully erosion, as well as their interactions over specific areas, particularly larger areas, underlining the importance of utilizing tools for scenario analyses. However, they also acknowledge that such tools generally have high data requirements, thereby complicating their application ability, particularly over larger areas. In these respects, the modelling framework proposed in this study, which benefits from the integration of satellite imagery and LANDPLANER, can be considered as an effective compromise between the need to have a representative approach for the complexity of gully processes and their dependence on main hillslope dynamics, and the need for easy application in different geo-environmetal contexts at meaningful scales.

| CONCLUSIONS
Based on the need for dynamically based predictions of the mechanism, timing, and occurrence of gully phenomena in high-erosivity landscapes affected by overcultivation, this study proposes a framework to analyse spatial and temporal gully occurrence, using freely available Sentinel-2 satellite data and the open-source software LANDPLANER. To improve over conventional methods which are limited by static geo-environmental input data, the proposed framework envisages the use of NDVI data for deriving dynamic CN maps used for the input of LANDPLANER, such that the model could predict gully occurrence using a twofold erosion modelling schema based, respectively, on erosion index and topographic threshold approaches.
The framework was evaluated and used to investigate the possible occurrence of gullies in different seasons for the Freddana Torrent catchment in the Tuscany region, Central Italy, for which a geomorphological gully inventory was available and used as a modelling benchmark. The framework proved to be effective to characterize the spatial and temporal occurrence of gully phenomena on a seasonal basis. In particular, the use of seasonal CN input maps derived from winter and autumn Sentinel-2 images performed better than the averaged CN input counterparts. The framework is applicable in all areas where topographical and LULC data are accessible, even if not highly detailed. Furthermore, the proposed methodology could easily be adapted to exploit very high-resolution commercial optical satellite images, for investigating the effects of specific rainfall events on gully erosion phenomena.
Although the method is slightly limited in terms of its dependence on the baseline gully inventory data, the overall outcome shows good promise for expanding to the application of multi-temporal and eventbased inventory maps for better characterization of small-scale dynamics.
The proposed method to convert remotely sensed NDVI data to CN values proved to be effective for providing a more reliable and realistic characterization of the conditions leading to runoff development and a better investigation of the spatial and temporal conditions triggering the gullies in the studied catchment. This method is expected to provide a scientific basis for more complex conversion schemas to be investigated in future studies.
The study was additionally funded by the Project BIOESSaNS were executed under the RStudio IDE, using the R language version 3.6.X.