Rainfall-Induced Landslides and Erosion Processes in the Road Network of the Ja é n Province (Southern Spain)

: Rainfall thresholds are one of the most widely applied methods for indirectly estimating landslide return periods, which are subsequently used in hazard analyses. In this study, the starting point is an incidence database of landslides and erosive processes affecting the road network of the province of Ja é n (southern Spain), in which the positions and dates of civil repair works can be found. Meanwhile, the use of a daily rainfall database in a dense grid (1 km) allowed for the estimation of the rainfall series at each incidence point with high precision. Considering the news in the local media and applying spatial proximity, temporal proximity, and maximum return period criteria, rainfall events of various duration (1 to 90 days) could be associated approximately with each point. Then, the rainfall thresholds and their return periods were estimated. A linear equation was adjusted for the rainfall duration threshold (E = 6.408 D + 74.829), and a power-law curve was adjusted for the intensity–duration pair (I = 47.961 D − 0.458 ). Non-signiﬁcant differences were observed between the thresholds and the return periods for the lower and higher magnitude incidences, but the durations for the former were lower (1–13 days), compared to those of the latter (7–22 days). From the equations, rainfall events of different durations could be estimated for use in hazard analysis, as well as for the future development of warning systems.


Introduction
Landslides are considered one of the most important natural hazards worldwide, causing thousands of victims per year and costs worth billions of euros [1][2][3][4][5][6]. Landslides originate in different mountainous regions in Europe [7], such as the Alps, Norway, and the Mediterranean countries [8,9], and specifically in Spain and the Betic Cordilleras [10][11][12]. Compared with other risk phenomena, such as earthquakes or floods, the effects of landslides are more diffuse and continuous in space and time; thus, according to some studies, their impact has been underestimated [4]. Despite this, they produce significant damage to infrastructure, properties, and the environment itself, as well as interrupting socioeconomic activity [1][2][3].
One of the most effective measures for risk prevention and mitigation is its evaluation which, according to the classic formulation of Varnes [1], includes both the hazard (probability that a potentially harmful phenomenon occurs in a given space and time) and the exposure and vulnerability of the elements at risk. For the former, there exist deterministic methods based on precise knowledge of the factors conditioning the instability [13]. These include geotechnical properties, terrain morphology, and the hydrological conditions usually conditioned by rainfall as the triggering factor [14][15][16][17][18]. However, the difficulty of obtaining accurate data [19], especially for studies carried out in a more or less extensive area, leads to probabilistic methods being more frequently applied in hazard studies. Probabilistic methods are based on correlation analyses between determinant factors and landslides, both by means of bivariate approaches [12,20,21] and multivariate ones [3,[21][22][23]; however, in recent years, machine learning techniques have become increasingly used [21,24]. The first step in these analyses is to elaborate inventories or databases that collect the spatial locations, occurrence times (dating), and thematic attributes of the movements [11,25]. For this, direct observation, geomatics data capture techniques, and in situ or laboratory tests of the different terrain properties are used [11,26].
Among them, dating is usually one of the most complicated issues for landslides, due to their aforementioned diffuse and continuous nature [27,28]. Direct dating requires recording by direct observation or in-situ sensors; however, geomatics techniques, such as Global Navigation and Satellite Systems (GNSS), photogrammetry, optical remote sensing (ORS), LiDAR, and Interferometry of Synthetic Aperture Radar (InSAR), have allowed for important advances, especially following their spatial and temporal resolution increases [29,30]. Another option is indirect dating from triggering factors [26,31], such as earthquakes and/or rainfall, which are more easily recorded by different instruments, usually gauges or digital sensors. Considering that rainfall is the triggering factor in most cases, it is necessary to establish the relationship between rainfall and landslides, which have been done in numerous studies worldwide [32][33][34][35][36][37][38].
In some works, very precise knowledge of the rainfall data (hourly), as well as the moment in which the movement starts, have been used [9,19,[41][42][43][44][45][46][47][48][49][50]53]; meanwhile, in others, only the daily rainfall data are known [28,36,40,51,52]. There are even cases in which the landslide time or date can only be approximated [36,40,43,51] and reconstructed after subsequent inventories and/or reviews from the news found in the media [28,40,51,52]. In any case, following [42], different variables can be used to define these thresholds, such as the total event rainfall (E), rainfall event-duration (E-D), intensity-duration (I-D), and rainfall event-intensity (E-I). Although most studies have been carried out at the local level, there are some cases of application over large regions [42][43][44][45][46], as well as works and databases that collect indices all over the world [43,54]. Moreover, some works have led to the development of algorithms and computational tools that can calculate these thresholds [48,49].
Once the thresholds are determined, they can be applied to predict the probability of landslide occurrence by determining the return periods [19,28,36,40,41,51,52], which can be incorporated into the corresponding hazard maps [64,65]. Likewise, they can be integrated into (early) warning systems [21,63,64], which can prevent the population and authorities from being subjected to landslides in those cases in which the rainfall threshold is reached.
The objective of this study is the determination of rainfall thresholds that cause landslides or erosion processes associated with the road network of the Jaén province. For this, an incidence inventory or database between 1997 and 2013 and a rainfall database between 1971 and 2016 are available. Then, news from the local media were also considered in the analysis that allowed associating incidences with rainfall events. This has led to the calculation of rainfall thresholds and the return periods, thus allowing for hazard modelling in the province of Jaén and the development of warning systems in future works.

Study Area
The study area (Figure 1) corresponds to the province of Jaén (13,486 km 2 ). It coincides approximately with the natural region of the eastern or upper Guadalquivir River basin. The altitude ranges between 152 and 2160 m, with an average value of 715 m. The average slope is 12.21 • , although it is also very variable between the mountain ranges and the lower lands of the Guadalquivir valley.  The predominant land-use is agricultural crops; within them, olive grove constitutes 44% (5928 km 2 ) of the province's surface area [69]. To a lesser extent, other crops (e.g., cereal) and areas of natural vegetation appear in the mountains, along with scrub and coniferous/hardwood forests.
The province has a population of 638,000 inhabitants, with only 2 urban areas From the geological and morphological point of view, three domains can be distinguished which, from South to North, are as follows (Figure 1b):

•
The External Zones of the Betic Cordillera, which are made up of mesozoic and cenozoic carbonate or loamy-clayey rocks, structured as a fold and thrust belt from the lower Miocene to the present [66], in which several paleogeographic domains appear (Prebetic and Subbetic). The Betic External Zone form several mountain ranges (Sierra Cazorla and Segura, Sierra Mágina and Sierra Sur of Jaén), partially isolated by main rivers and tributaries of the Guadalquivir River.

•
The sedimentary infill of the Guadalquivir basin, differentiated into two parts. In the north, the Guadalquivir basin is filled with Miocene loamy and clayey sediments, which are slightly deformed and which overlie the tabular cover of the Iberian Massif, made up of Triassic clays and sandstones and Jurassic limestones. In the south, the infill of the basin is highly deformed by the Betic Miocene displacements, which incorporate tectonically Betic soft materials as Triassic evaporites (salt and gypsum) or Cretaceous clayey marls [66].

•
The Variscan Domain, which constitutes the outcropping basement of the Iberian Massif, in which metapelites (slates, grauwackes, and so on) and intruding igneous rocks (granites and granodiorites) are the predominant lithologies.
Over all these materials, quaternary deposits related to present fluvial dynamics and slope sediments are located.
From the climatic point of view, the province of Jaén corresponds mostly to the hot summer Mediterranean (Csa de Koppen) [67] climate type. More specifically, it can be catalogued as the Mediterranean meridional type of the Guadalquivir valley [68]. This is characterized by a mean annual precipitation (MAP) between 500 and 650 mm, with maximum values distributed between the autumn, winter, and spring and minimum values in the summer. Meanwhile, the average temperatures are 17-18.5 • C, with very pronounced maximum values in summer. However, there are sectors in the province with MAP higher than 1000 mm in the mountain ranges (Mediterranean mountain), and others with MAP that does not reach 400 mm in the southeast (Mediterranean arid), as shown in Figure 1c.
Thus, at a central point representative of the average physical conditions of the province (point 066 of the incidence database, see below), the MAP was 533 mm within the period considered (1971-2016), with a minimum value of 223 mm in the hydrological year 2004-2005 and a maximum value of 1026 mm, which was reached in 2009-2010 ( Figure 1d). This wide interval shows the variability of precipitation over the years, with a standard deviation of 172 mm and coefficient of variation of 0.32. Within the year, rainfall was higher between November and April (50-60 mm) and lower between June and September (below 25 mm), with a monthly average rainfall of 41.4 mm (Figure 1e).
The predominant land-use is agricultural crops; within them, olive grove constitutes 44% (5928 km 2 ) of the province's surface area [69]. To a lesser extent, other crops (e.g., cereal) and areas of natural vegetation appear in the mountains, along with scrub and coniferous/hardwood forests.
The province has a population of 638,000 inhabitants, with only 2 urban areas exceeding 50,000 inhabitants [70], and a low industrial activity focused on agriculture. On the other hand, it has a road network of different orders (state, regional, and provincial), in which the A-44 and A-316 highways stand out. This study is focused on the extensive and penetrative road network of the Provincial Deputy of Jaén, which is about 1600 km long. This network is a fairly representative sample of the different physical environments of the province and, so, its study can provide valuable information on the instability conditions, not only in the network environment, but also in the whole province.

Incidence Database
The overall methodology followed in this study is shown in Figure 2. The first step is the elaboration of the database of incidences on the road network of the Jaén province. This includes data extracted from the files of the works carried out for their repair and maintenance, completed with field data and other data extracted from previous maps. The database has been elaborated through the years, by gathering information on road interventions that took place from 1998 to 2013.

Rainfall Data Processing
Different meteorological databases were used to estimate the rainfall series: • Spain02 Database, high-resolution daily precipitation data, developed by the Institute of Physics of Cantabria (Spain) and the Spanish Meteorological Agency (AEMET) from a dense network of more than 2500 quality-controlled stations for precipitation and near 250 for temperatures. The Spain02.v5 provides daily data from 1951 to 2015, gridded in increments of 0.1°, corresponding approximately to a resolution of 10 km [71,72]. • RIA database, a network of agroclimate information by the Department of Agriculture, Fisheries, and Rural Development of the Andalusian Government [73]. It contains updated data on the networks of automatic meteorological stations (~120 stations), which are equipped with electronic sensors and distributed throughout the Andalusian territory. The data processing comprised the integration of the previous databases and the application of physical and statistical filters, which allowed us to obtain an interpolated regular grid at 1 km. This grid was used to assign a daily rainfall value, from 1971 to 2016, to each of the 186 incidence points of the road network in the province of Jaén. The value assigned was that of the closest grid node which, taking into account its resolution, is a Thus, the original database included data identifying the incidence (coordinates, project, works, and so on); vegetation and land-use; geomorphological and topographic data; geotechnical data (load capacity, constructive conditions, and so on); hydrogeology (drainage and permeability); description of the incidence (year, month, typology of incidence, roads, kilometres, and so on); geology (lithology and surface formations); road data (road surface, slope, curvature radius, and so on); and, finally, the constructive solution adopted. The database was subsequently tested on the ground, especially with regard to geology and the descriptions of the incidences, as well as morphological aspects.
Finally, the recorded and reviewed incidences were digitized onto the orthophotography and subsequently refined using several GIS tools (e.g., clipping, buffering). The result was an enriched incidence inventory or database of the road network of the province. A basic distinction was made between two categories: the lower magnitude and shallower processes that affect the road surface, road cuts, and embankments; and the higher magnitude processes that involve a certain general slope instability. Among the former, the following types were differentiated: erosive processes (gullies), undercut of road embankments, and small slides and collapses of the road cuts. Among the latter, slides, earth or mud flows, and creeping processes were differentiated. The inventory is presented and described in Section 3.1.

Rainfall Data Processing
Different meteorological databases were used to estimate the rainfall series:  [75].
The data processing comprised the integration of the previous databases and the application of physical and statistical filters, which allowed us to obtain an interpolated regular grid at 1 km. This grid was used to assign a daily rainfall value, from 1971 to 2016, to each of the 186 incidence points of the road network in the province of Jaén. The value assigned was that of the closest grid node which, taking into account its resolution, is a point located at a small distance (lower than 1 km). For the daily precipitation values in the province of Jáen, three zones can be considered, in two of which the WRF underestimates the precipitation and, in the other, it is weakly overestimated. To reduce the uncertainty of the data provided by the WRF, we used data obtained at ground stations as control points. Using geostatistical techniques, the uncertainty was less than 12% in all cases. Therefore, although the densification of the grid tended to smooth the real values, the estimation error and the uncertainty derived was low.

Rainfall Event Identification
The identification of rainy events associated with the landslides and erosive processes in the study area was based on relating the incidence database with the rainfall series.
First, from the daily data, the accumulated rainfall over 2 days, 3 days, 5 days, 7 days (1 week), 10 days, 15 days, 30 days (1 month), 45 days, 60 days (2 months), 75 days, and 90 days (3 months) were calculated, in order to analyse the influence of short-and medium-term rainfall on the generation of landslide and erosion processes.
Then, different rainfall variables could be defined: rainfall amount associated with the event (E) in mm; duration of the event (D) in days; and the intensity (I), which was calculated as the relationship between the rainfall and the duration of the event and expressed in mm/day. For each rainfall event-duration (E-D) pair, the probability of exceedance and return period (T) in years were calculated considering a Weibull series, as in previous studies [36,40,51,52].
Meanwhile, the incidence database only included information about the month in which the civil work started, while the accurate date when the incidence occurred remained unknown. Thus, additional information was used, in order to estimate a more precise date. For this, we used information found in the media, especially the news published in the local and regional press, such as the IDEAL newspaper [76], which has a historical record since 2006 and can be accessed freely on the internet. This approach has been used in some previous works [28,40]. After a deep search, based on terms such as landslides (and their synonymous terms in Spanish), road affected, traffic interruption, and so on, a total of 98 news items were found between 2006 and 2013, of which 27 events were directly related to the incidences (i.e., landslides or erosive processes). These are summarized in Table 1. Then, three criteria were applied: the spatial proximity, the temporal proximity, and the magnitude of the rainy event.

1.
The spatial proximity between the approximate location in the media and the incidence coordinates were estimated in the GIS. First, sections of roads and affected towns or municipalities mentioned in the news were selected. Then, a spatial query allowed for the identification of those incidences close to them. Five classes were established, depending on the distance: Class 1, 0-1 km; Class 2, 1-2 km; Class 3, 2-5 km; Class 4, 5-10 km; and Class 5, more than 10 km. 2.
The temporal proximity was addressed by analysing the time interval between the date of the news appearing in the local media and the month associated with the previously selected incidences. Five classes were also considered: Class 1, 0-3 months; Class 2, 3-6 months, Class 3, 6-12 months, Class 4, 12-24 months; and Class 5, more than 24 months. Summing the classes for the spatial and temporal proximities, only those incidences with a maximum of 6 points (e.g., Class 3 in both, or Classes 2 and 4 in each one) were selected. Thus, each incidence point could be associated with several rainfall events and their rainfall-duration (E-D) pairs. In addition to the events identified from the news, the complete rainfall series for the two years previous to the month of each incidence were examined, searching for the major events in each interval of duration. If events different from the above were found, they were also added to the database.

3.
Finally, the magnitudes of the rainfall events were considered. First, following some previous studies [36,40,51], the E-D pair with the longest return period was selected, for each of the events associated with an incidence, as the one most likely to trigger it. Moreover, according to [40], of all the possible events (E-D pairs) associated with each incidence, those which presented a return period of fewer than five years were discarded, as they were considered non-relevant for incidence triggering. This procedure allowed for enrichment of the incidence database, thus including several E-D pairs for each incidence (see some examples in the results Section 3.2).

Rainfall Threshold Calculation
Prior to the determination of the thresholds, the rainfall variables were analysed, by calculating their mean and modal values, both globally and for each of the typologies considered, in order to determine whether there were differences between the landslides of different typologies and magnitudes.
In this work, the determination of thresholds of the rainfall-duration (E-D) type was considered, which usually respond to linear equations of the type: although equations with a power-law can also fit: This type of threshold has been considered more appropriate for cases where only daily data are available [28,51,52]. They allow for knowledge of the amount of rainfall necessary to generate landslide or erosive processes, depending on the number of days. Nevertheless, thresholds of the intensity-duration (I-D) type were also determined, although these are more commonly used when intensity per hour (mm/hour) data are available [42,43]. In this case, power-law equations were adjusted. Both thresholds were calculated globally for all the incidence points, but they were also discriminated by typologies.

Incidence Database
The incidence database is shown in the map of Figure 3. It shows the typology and magnitude of the incidences, according to published classifications of landslides [77,78], and includes some significant examples in the study area. In general, practically all the incidences corresponded to shallow phenomena but, within them, two types were differentiated, depending on their size or magnitude [79]: • Very shallow processes, with magnitude between extremely and very small (<5000 m 3 ). These correspond to ruptures in the road cut, either of the slide or collapse-rockfall typologies, but also undercuts of the road embankment. Meanwhile, erosive processes (gullies) were identified, which also produce incidences on the roads. • Shallow processes in which there is mobilization of the slope where the road is located, with a magnitude generally between small and medium (5000-500,000 m 3 ). Within these, slope movements of a slide or flow type were considered, according to [77,78]. Soil creeping processes were also distinguished from those flows which were welldefined in the landscape.
The distribution by typology is shown in Table 2. As can be observed, there were 46 incidences corresponding to gully processes, 47 punctual incidences in road cuts (38 landslides and 9 collapses), and 30 incidences associated to undercutting in road embankments. Then, 77 incidences were directly related to the road; that is, to human activity. Moreover, there were 63 landslides of higher magnitude, among which 21 slides, 26 flows, and 16 creeping areas were identified.   Meanwhile, Table 3 shows the distribution of the year in which the civil works to repair the road started. Most of them were concentrated into two years, 2010 (61 incidences) and 2013 (70 incidences). Table 3. Distribution of incidences by year in which the civil work started.

Year
Number Year Number   1998  2  2006  2  1999  5  2007  1  2000  3  2008  2  2001  3  2009  3  2002  2  2010  61  2003  1  2011  18  2004  3  2012  3  2005  4  2013  70 3.2. Rainfall Events Figure 4 shows the daily rainfall series associated with the aforementioned significant incidence points (shown in the map of Figure 3), where different rainfall events can be observed. Following the methodology described, the rainfall events associated with each incidence were searched in the two years (24 months) previous to the starting of the repair work. Thus, Figure 5 shows the two-year rainfall series for different event durations in two significant incidence points (099 and 181, not shown in Figure 3). Some arrows in red have been included to point out the rainfall events identified in each incidence, which were later used in the thresholds calculation.
A total of 446 rainfall events (E-D pairs) associated with the 186 incidence points were found that met the established criteria of spatial and temporal proximity, as well as the maximum return periods. Thus, 17 points were associated with 4 potential events, 60 points with 3, 86 points with 2, and 21 points with a single event. Meanwhile, some events affected the whole province in a general way and, therefore, the provincial road network, while others affected more restricted sectors. Some of them, those affecting a minimum of 5 points, are shown in Table 4.     The second important rainy period occurred in the year 2012-2013, with several events: 03-08/11/2012, with daily rainfall that exceeded 50 mm, 2-day rainfall around 70 mm, and weekly rainfall that reached 150 mm (Figures 5b-d and 6e,f), all generalized in the road network (50-60 points). The return period was 9 to 15 years. Subsequently, on 18/03/2013, rainfall of 166 mm was reached in 15 days at 43 points, with a return period of 6 years.
More locally, there were other rainfall events potentially associated with incidences such as those at the end of 1997, when 400 mm was exceeded in 60 days (Figure 6a Table 5 shows the mean values of the rainfall amount, duration, and intensity (E, D and I) variables, as well as the modal value of the duration, the total events, and discriminated by typology. Likewise, Figure 7 shows the histograms of the duration of the events. From the data shown in Table 5 and Figure 7, a shorter duration of rain events was generally observed in the lower magnitude incidences, with an average value of 13.70 days (the modal value being 1 day). Meanwhile, for the higher magnitude incidences, such as slides, flows, and creeping processes, the average duration was 21.60 days with a modal value of 7 days. In the lower magnitude incidences, the average rainfall was 223 mm and the intensity was 22 mm/day. For the higher magnitude landslides, the average rainfall was 155 mm, with an average intensity of 28 mm/day.

Rainfall Thresholds
Regarding the thresholds, Table 6 shows the equations obtained both for the rainfallduration (E-D) threshold (linear and power-law adjustment), and for the intensity-duration (I-D) threshold (power-law adjustment). The table also shows the coefficient of determination (R 2 ) of the adjustment. Figure 8 shows these thresholds for lower and higher magnitude incidences.
As can be seen from Table 6 and Figure 8, the equations were quite similar for both the lower and higher magnitude incidences; although, in the case of linear adjustments, the intercept in the former (72)-especially in the collapses (61)-was lower than that in the latter (82), the slopes being similar. In the same way, the power-law base was somewhat lower in the lower magnitude incidences (47.5) than in the higher magnitude ones (49.7), both for the E-D and I-D thresholds. Table 6. Equations for E-D (linear and power-law) and I-D thresholds.

Discussion
From the results obtained, it can be observed that the road network in the Jaén province was affected by numerous incidences in the period studied. A set of 186 incidences was registered, of which 123 (66%) corresponded to lower magnitude processes (gullies, undercut of road embankments, slides and collapses in road cuts,) and 63 (33%) corresponded to processes that were also shallow but of higher magnitude, affecting the entire slope (slides, flows, and creeping). Although all the movements affected the road network, some of them (undercut on road embankment, slides and collapses in road cuts, comprising 47% of the incidences) can be considered as directly related to human activities.
The temporal distribution was quite irregular, concentrated mainly in two years: 2010 (with 60 incidences) and 2013 (with 71), representing more than 75% of incidences.
This distribution seems to be related to the occurrence of rainy periods in the province. Without the contribution of other factors, given the low tectonic activity [80], rainfall was the main triggering factor for landslides in the province. As indicated in the introduction, the relationship between landslides and rainfall has been well-established throughout the world [32,33,35], particularly in Europe [42,43] and in the Mediterranean countries [8,9,44,46]. In Spain, these relationships have also been found [28,51,52,81].
Thus, the simple observation of the distribution of the mean annual precipitation (MAP) of point 066, representative of the whole province (Figure 1d), as well as the series of daily rainfall associated with different incidence points (Figure 4), allowed for the establishment of this relationship. Thus, considering the mean annual precipitation (MAP) of 533 mm at point 066, 1026 mm was reached in the hydrological year 2009-2010, 911 mm in 2010-2011, and 950 mm in 2012-2013. The activity of these years has also been observed in natural slopes of some sectors of the province [30]. Meanwhile, the remaining years barely exceeded 600 mm, except for 1976-1977, 1995-1998, and 2003-2004. However, discarding the first years in which there was no recording of incidences, the intense rainfall of other years, such as 1997-1998 and 2003-2004, was not reflected in the incidence database, as explained by the lower magnitude of the rainfall events or because the incidences were not registered (being in the first years of the database elaboration).
Analysis of the rainfall series associated with each incidence made it possible to more precisely identify a set of possible events for each point in different intervals of antecedent rainfall (duration), based on criteria of spatial and temporal proximity, in relation to the news that had appeared in the local media. Figure 5 shows some of these possible events with different durations, associated with two significant incidence points of different magnitude. Thus, for the point 099 (a higher magnitude incidence), several events were observed in the hydrological years 2009-2010 and 2010-2011; while, for the point 174 (a lower magnitude incidence), several events took place in the hydrological year 2012-2013. Some of these events corresponded to those that occurred with different magnitude in some sectors of the province, or those affecting the whole of the road network in a generalized way, as shown in Table 7. They also coincided with the years in which the MAP was higher, as mentioned above.  At the point 099 (whose civil works started in June 2011), the two-year antecedent rainfall series began in the hydrological year 2009-2010, with a rainfall of 150 mm in 7 days to more than 200 mm in 15 days in December. Rainfall continued in January, reaching more than 300 mm in 30 days and, in February-March, it exceeded 570 mm in 75 days and 680 mm in 90 days, more than the annual precipitation in most of the province. After summer 2010, the rainfall recovered, reaching values of 120 mm in 7 days, 230 mm in 15 days, and near 340 m in 30 days for December. Thus, four E-D pairs were selected, for different dates that met the aforementioned criteria ( Figure 5), with return periods higher than 20 years (even reaching 45 years), as the most important events of the entire rainfall series. Meanwhile, at point 174 (whose civil works started in June 2013), the series of antecedent rainfall reached important values from September and, especially, November 2012, with events of daily rainfall greater than 60 mm and weekly rainfall of 150 mm, accumulating about 300 mm in 30 days by the end of this month. After this rainy period, the rainfall decreased, but recovered in March 2013 and reached values close to 130 mm in 7 days. Then, four E-D pairs were selected, with return periods always greater than 5 years (mostly between 10 and 20 years). These two particular examples coincided with the following analysis, in which the lower magnitude incidences were usually related to intense rainfall of short duration (1-7 days), while higher magnitude incidences required a longer duration (1-3 months).
Analysis of the average values of the considered variables showed some aspects of interest, such as a shorter duration of the events associated with lower magnitude incidences (mode of 1 day and mean of about 13 days), compared to those of higher magnitude (mode of 7 days and mean higher than 20 days). Consequently, the amount of rainfall was lower in the former (around 150 mm) than in the latter (around 225 mm), unlike the intensity (28 and 22 mm/day, respectively). This difference in behaviour has been pointed out in previous studies considering deep landslides [36,40], where the duration of antecedent rainfall ranged between less than 15 days for shallow landslides to more than 30 days for deeper ones [36,40]. Usually, the prediction is more complex for deeper landslides, for which it is necessary to consider the antecedent rainfall that determines the soil moisture conditions in the medium-term [16,18,19,28,[40][41][42][43]51,52,57,[59][60][61][62]81], or even the variation in annual rainfall over several years [58]. These landslides are triggered by a reduction in the shear strength of affected soils and rocks, related to the constant increase in groundwater level as an effect of long-term rainfall periods [40]. Although deep landslides were not considered in the strict sense in this study, a certain difference was observed between the lower and higher magnitude incidences, with respect to the duration of events causing the incidences.
In any case, the consideration of antecedent rainfall provides a simple way to introduce hydrogeological conditions into these studies, even when only shallow landslides are analysed. Some studies, which have mainly used daily data (such as in this one), have been based on the analysis of rainfall duration periods longer than 1 month [51], while other, more sophisticated methods distinguish between the antecedent and triggering rainfall [28,52]. The use of a calibrated antecedent rainfall that decreases with time [38,40] provides another way to simulate the hydrological conditions; however, in this preliminary study, in which the incidence date had high uncertainty, a simple method based only on rainfall amount-duration pairs was applied.
Regarding the rainfall thresholds, we focused especially on E-D type pairs, as they are two truly independent variables [45,46,56]. In addition, in this case (as in other ones), we utilised daily data [28,40,51,52,62,82], which can make it difficult to accurately calculate the intensity, with respect to those other studies that use hourly data [42,43]. Nevertheless, the adjustment of linear equations was tested, with good results, as in previous studies with daily data [51], as well as power-law curves for both the E-D and I-D thresholds ( Table 6 and Figure 8). The different thresholds show a good general fit (R 2 higher than 0.9), without significant differences between them. Moreover, normalized thresholds using the MAP were calculated, and all the equations were expressed using the duration (D) in hours in order to facilitate comparison with other studies. Most of these results were not included in the previous section, for simplicity, but are summarized in Table 7, where a comparison with the thresholds obtained by other authors is shown.
In general, the thresholds were within an order of magnitude of those corresponding to other studies, validating the results of this work. However, it must be taken into account that some of them used lower value thresholds; that is, thresholds adjusted such that only a reduced percentile (10, 5, or 1%) of the landslides is under the threshold curve. In this study, we mainly considered thresholds adjusted to the mean values, but the threshold of lower values (5%) was also calculated for the linear adjustment (Table 7, low). From both thresholds, the rainfall amount and intensity corresponding to the different durations were calculated, which could be used to develop a warning system regarding the activation of incidences in the road network of the province of Jaén (Table 8). Hydro-meteorological thresholds [59][60][61][62] are increasingly being used to overcome the drawbacks of rainfall thresholds, which do not consider adequately the hydrological conditions of slopes and, besides, produce a great proportion of false positives in the prediction, thus limiting the development of warning systems. There exist different methods to determine these thresholds, some of them based on hydro(geo)logical models which take into account detailed data of the phreatic level, soil humidity, porosity, permeability, saturation, rainfall, and so on [59][60][61], while others are based on data at the basin level, such as rainfall, evapotranspiration, runoff, and so on [62]. In both, additional data are necessary, which were not available in this case, such that they were discarded herein.
Considering the differences by typology, only small differences were observed between the values obtained; the intercept of linear adjustment and the base of the power-law being somewhat higher in the incidences of higher magnitude than in the lower magnitude ones. It must be taken into account that, in any case, the incidences affecting the road network were always shallow, and moreover, the temporal resolution of the data (daily rainfall) was most likely too low to find differences. Differences between the incidences more directly related to the human activities regarding to those less related were not observed. The only typology that presented a certain difference with respect to the remaining incidences, were the collapses (with a lower slope and intercept). This typology is usually associated with steeper areas and road cuts (which are prone to landslides) and the thresholds are likely lower than in other incidences. However, as discussed above, the duration showed higher values in the higher magnitude incidences then in the lower magnitude ones.
The importance of assessing the uncertainty in the determination of rainfall thresh-olds has been discussed in several studies, given that they are based on empirical data which are not always acquired with the required accuracy [42][43][44], both in the spatiotemporal component and in the measurement component (rainfall gauges). Approaches based on bootstrapping techniques [45] have been proposed, in order to estimate the uncertainty of the power-law parameters (α as the base or scaling parameter; and β as the exponent or the shape parameter) in the E-D thresholds. These approaches have been used in other studies oriented to the development of warning systems [45,46,48,49,53], where uncertainty values representing 5-10% of these parameters have been found [50]. The influence of the temporal resolution on the uncertainty has also been analysed, resulting in smaller scaling parameters (intercepts in logarithmic scale), higher shape parameter (slopes), shorter ranges of validity of the thresholds, and higher uncertainties when the temporal resolution decreased [53]. In this study, the main uncertainty source was not the spatial component, as the meteorological data grid was quite dense (1 km) but, instead, the temporal component (daily rainfall resolution), and especially the uncertainty of the incidence date. This date was estimated from the month when the civil work started, applying the criteria of spatial and temporal proximity to the news appearing in the local media, as well as the magnitude of the event (return period). In these conditions, a validation analysis with the use of positive and negative events is difficult to implement. These limitations, in terms of the data, and the preliminary nature of the study, led us to not consider the uncertainty analysis. Nevertheless, future works utilising more data and with the objective of developing a warning system should address these concerns regarding validation and uncertainty analyses.
Finally, a brief note on the relationships between landslide activity, rainfall, and the global climate could be necessary, although it exceeds the objectives of this study. Several studies [36][37][38]40,52] have pointed to a relationship between rainfall of high intensity (which generates landslides) and global climatic phenomena. Among these, the wellknown teleconnections stand out, such as the South Pacific-El Niño oscillation (ENSO), the North Atlantic oscillation (NAO), or the Western Mediterranean oscillation (WeMo). All of them are related, and may have an influence on the rainfall of the Iberian Peninsula and the Mediterranean [84,85], as has been observed in the Balearic Islands [37].
It is well-known that the high negative anomalies of the NAO (NAOi) are the origin of rainy winters and storms in Portugal and the whole of the Iberian Peninsula [84], thus producing floods [86] and landslides [28,36,40,52]. In this case, they are mainly related to deep landslides [40], but also to shallow landslides [36,52]. Meanwhile, the WeMo also seems to influence the rainfall in the Mediterranean area of the Iberian Peninsula [87]. Thus, in southern Spain, the relationships between both indices and rainfall have been analysed [52]. Although a significant correlation has not been obtained, two of the wettest hydrological years, such as 1995-1996 and 2009-2010 (the latter being one of the years with the highest activity in the study area), were observed to be related to negative anomalies of both indices.
Regarding the role of climate change in rainfall-triggered landslides, this is an interesting issue currently in discussion [88], although the prediction of rainfall events that generate landslides has shown a high level of uncertainty [89]. However, it seems clear that a global climate change scenario, in which severe events such as intense rainfalls are expected [90], should influence landslide activity [91], specifically the rainfall thresholds [88]. Predictions for Spain and the Mediterranean [92,93] have pointed out, in this direction, that Severe Weather Threats (SWEATs) are estimated to increase in the near future, especially in the summer and autumn [94]. Thus, higher levels of hazard and risk affecting civil infrastructure in the province of Jaén are also expected, as has been shown by a preliminary approach, in which a 30% increase in daily rainfall was found in some sectors of the province [95].

Conclusions
Rainfall thresholds are one of the most widely applied methods for indirectly estimating landslide return periods, which are subsequently used in hazard analyses. They are also used in the implementation of early-warning systems.
In this work, the starting point was an incidence database in the road network of the province of Jaén, in which the positions and dates of civil repair works are contained. Moreover, a meteorological database with daily rainfall data in a dense grid (1 km) allowed us to link the rainfall series to the incidence points accurately. Then, the identification of rainfall events that potentially generated landslides and erosion processes was addressed, using criteria related to the spatial proximity, temporal proximity, and return period, additionally considering the news appearing in the local media.
Regarding the results, some relevant aspects could be extracted: • Several events were identified, the most important being related to the hydrological years 2009-2010 and 2012-2013. Some of them were located in specific areas and other ones affected practically the entire road network. The return periods of these significant events were always greater than 5 years and, in some cases, exceeded 10-20 years.

•
The lower magnitude incidences usually presented a shorter duration (mode of 1-15 days), compared to those of higher magnitude (7-30 days). Consequently, the amount of rain was lower in the former (around 150 mm) than in the latter (around 225 mm).

•
The thresholds obtained for both the rainfall-duration (E-D) and intensity-duration (I-D) pairs were on the same order of magnitude as those calculated by other authors, some of them in a similar environment (i.e., Mediterranean countries). The different types of thresholds tested (E-D or I-D, linear or power-law) showed a good fit, without significant differences, likely due to duration data being in units of days, not in hours, and the shallow nature of all the incidences.

•
In this case, there were no differences in the thresholds between the lower and higher magnitude incidences, unlike the variables (E, I, D) themselves. • Finally, from the thresholds, rainfall amounts and intensities for different durations of the events were calculated (e.g., about 80 mm for 1 day and more than 250 mm for 1 month), considering not only the threshold adjusted to the mean values but the threshold adjusted to the lower values.
Future improvements to the study should first address extension of the database to other roads, and even the inclusion of landslides on natural slopes, which would allow the thresholds and calculated variables to be refined. In this sense, recording the incidences or landslides with precise knowledge of the date of occurrence could contribute to this refinement. The use of more accurate data would also allow for addressing advanced models, in which antecedent and triggering rainfalls are distinguished and, even (if additional data were available), the determination of hydro-meteorological thresholds. Another future approach consists of estimating the uncertainty of the models and their validation (e.g., using different samples of incidences for training and testing, or by means of temporal validation). Then, these thresholds may be reliably used in hazard analysis or for the implementation of (early) warning systems in the study area.