Deformation responses of slow moving landslides to seasonal rainfall in the Northern Apennines, measured by InSAR

Slow moving landslides are widespread geomorphological features in the Northern Apennines of Italy where they represent one of the main landscape forming processes. The lithology of the Northern Apennines fold and thrust belt is characterized by alternations of sandstone, siltstone and clayshales, also known as ﬂysch, and clay shales with a chaotic block in matrix fabric, which are often interpreted as tectonic or sedimentary m´elanges. While ﬂysch rocks with a high pelitic fraction host earthslides that occasionally evolve into ﬂow like movements, earthﬂows are the dominant landslide type in chaotic clay shales. In the present work, we document the kinematic response to rainfallof landslides in these diﬀerent lithologies using radar interferometry. The study area includes three river catchments in the Northern Apennines. Here, the Mediterranean climate is characterized by two wet seasons during autumn and spring respectively, separated by dry summers and winters with moderate precipitation. We use SAR imagery from the X-band satellite COSMO SkyMed and from the C-band satellite Sentinel 1 to retrieve spatial displacement measurements between 2009 and 2016 for 25 landslides in our area of interest. We also document detailed temporal and spatial deformation signals for eight representative landslides, although the InSAR derived deformation signal is only well constrained by our dataset during the years 2013 and 2015. In spring 2013, long enduring rainfalls struck the study area and numerous landslide reactivations were documented by the regional authorities. During 2013, we measured higher displacement rates on the landslides in pelitic ﬂysch formations compared to the earthﬂows in the clay shales. Slower mean velocities were measured on most landslides during 2015. We analyse the temporal deformation signal of our eight representative landslides and compare the temporal response to precipitation. We show that earthslides in pelitic


Introduction
In the Northern Apennines, slow moving landslides are among the most common landscape forming agents, in particular where mechanically weak rocks form the bedrock (Bertolini et al., 2004;Simoni et al., 2013).According to Bertolini et al. (2004), most of these landslides started to form during the early holocene and continue to exhibit slow deformation.However, this slow motion may be interrupted by periods of high displacement rates, which is why these landslides represent a major natural hazard, that periodically causes damage to houses and infrastructure (Bertolini et al., 2005).Synthetic aperture radar interferometry (InSAR in the following) provides the possibility to assess deformations of slow moving landslides both on different scales and the technique was repeatedly used to infer spatio-temporal landslide displacement in different parts of the world (Colesanti et al., 2003;Hilley et al., 2004;Roering et al., 2005;Iasio et al., 2012;Handwerger et al., 2013;Wasowski and Bovenga, 2014;Handwerger et al., 2015).Handwerger et al. (2013Handwerger et al. ( , 2015) ) successfully measured seasonal deformation of large earthflows in a Northern California river catchment, that is characterized by mélange type shales similar to those in the Northern Apennines.In this paper we present the results of an InSAR campaign that contain datasets derived from the X-band satellite constellation COSMO SkyMed (CSK in the following) and the C-band satellite Sentinel 1A (Se1 in the following).CSK has a spatial resolution of ca 3 m and a variable acquisition frequency, that in our case range from 16 to 30 days, while Se1 has a moderate spatial resolution of ca.20 m and a high acquisition frequency of ca. 12 days.In this paper, we document the spatial deformation patterns of 25 landslides during the period 2013-2016 in three river catchments covering an area of ca.3600 km 2 .Flysch rocks with high fraction of pelitic material (pelitic flysch in the following) and chaotic clay shales, both of which belong to the Ligurian Nappe of the Northern Apennines (Bettelli and Panini, 1992;Pini, 1999), are the two main lithologies that are affected by landsliding (Bertolini et al., 2005).During the spring 2013 the study area was struck by long enduring rainfall that caused numerous landslide reactivations (see also Pizziolo et al., 2013Pizziolo et al., , 2015)), whose displacement exceeded tens of meters in some cases (Berti et al., 2017).Overall, 92 landslides were registered by the regional authorities during spring 2013.79 landslide reactivations were registered in pelitic flysch formations while 9 occurred in chaotic clay shales and only 4 were counted in other lithologies like coarse grained flysch, pure sandstones or marls (see Fig. 1 in the supplementary material for a more detailed representation).
The present study was designed in order to test if this different behaviour of the landslides in the two main lithologies is also observable in the InSAR derived displacement signals.We analyse morphological differences and similarities between landslides that yielded InSAR signals in different bedrock lithologies and inspect the relationship between their displacement history and meteorological forcing.Seasonal trends and reactivation episodes are documented and discussed in detail for eight selected cases.We show that landslides in pelitic flysch formations exhibit temporal deformation patterns that correlate with periods of long persistent rainfall, while only weak seasonal trends are detected for landslides in clay shales.

Geological and geographical background
The Northern Apennines are a fold and thrust belt in Italy that formed due to the convergence of the European and Adriatic plates (Boccaletti et al., 1971;Picotti and Pazzaglia, 2008).In the study area, which is located in the Northern Apennines South of Bologna and Modena, chaotic clay-shales with block in matrix fabric (Pini, 1999;Vannucchi et al., 2003) and deposits of turbidity currents, also known as flysch (Ricci Lucchi, 1986), are the most common lithologies (Fig. 1 a).In particular the tectonically shearded formations of the Ligurian Unit have a high fraction of pelitic material, which is why we refer to them as pelitic flysch.In zones where chaotic clay shales form the bedrock, earthflows with distinct source, transport and deposition zones are the dominant type of slope failure (Simoni et al., 2013).Because the pelitic flysch formations and the chaotic clay shales have poor mechanical characteristics they represent the lithologies in our study area that are most affected by landsliding (Fig. 1,b, c and Fig. 1 of the supplementary material).Where the slopes are composed of pelitic flysch, a broader spectrum of landslide types can be observed, but often the failures have characteristics of translational or rotational earth-or rockslides in the upper portions of the landslide and may propagate as earthflow further down-slope (Borgatti et al., 2006;Corsini et al., 2006;Berti et al., 2017).Bertolini et al. (2004) carried out radiocarbon dating on tree logs, stumps and deposits of peat bogs that were buried inside landslide deposits of the Northern Apennines during catastrophic failures.They found ages ranging from the transition between Pleistocene and Holocene (ca.13.759 years B.P.) and Medieval times (ca.500 years B.P.), with most landslide ages grouping after the transition between the Atlantik and the Sub boreal periods (4900-3900 years B.P.).
Today the Northern Apennines have a Mediterranean climate and total annual precipitation reaches on average 1300 to 1400 mm (Berti et al., 2012).The temporal pattern is characterized by intense rainfall in spring and autumn, separated by dry summers and winter months with moderate precipitation that in part occurs as snowfall (Tomozeiu et al., 2000(Tomozeiu et al., , 2002;;Pavan et al., 2008) However, in our study area total spring rainfall accumulated to 532 mm in 2013, 426 mm during 2014 and 384 mm during 2015.Also maximum daily intensities were with 41 mm/day (2013), 40 mm/day (2014) and 39.181 mm/day (2015) higher than those during 2000-2016.Fig. 2 in the supplementary material illustrates that in our study area more rainfall was measured in greater heights and that little variation occurs for areas of equal height.Berti et al. (2012) used bayesian inference to develop a practical rainfall threshold for landslide triggering in the study area, taking into account both rainfall events that caused strong reactivations as well as events that did not.They showed that landslide probability increased with higher rainfall duration and intensity, showing no distinct limit between critical and not critical rainfall.
While antecedent rainfall appeared not important for landslide triggering itself (Berti et al., 2012), Fig. 1: a) Lithological map of the study area (Berti et al., 2012) with the reactivations that were registered by the regional authorities during spring 2013 and the raingauges that are available in the study area.b) Lithological map with landslide perimeters over the western part of the study area, where pelitic flysch formations crop out more extensively compared to c) the eastern part of the study area, where chaotic clay shales are more common.The landslide perimeters area derived from a detailed regional inventory (Servizio Geologico, Sismico e dei Suoli della Regione Emilia-Romagna, 2014) and In Fig. 1 of the supplementary I we report the main statistical aspects of this landslide inventory, the reactivations and the annual rainfall derived from the raingauges pattern in the study area.
the relationship between deformation and rainfall is often only clear when longer periods of preceding rainfall are taken into account (Borgatti et al., 2006;Ronchetti et al., 2007).Moreover surface deformation varies often for different areas of the landslide (Iverson and Major, 1987;Petley et al., 2005) and the hydrological, hydrogeological conditions on a slope can be complex (Cervi et al., 2012).Anisotropy in primary conductivity, but also the presence and the formation of fissures and macropores on a landslide are difficult to asses on a slope scale, but may influence significantly infiltration, water flux, drainage and hence the build-up of pore pressure in the landslide body (Krzeminska et al., 2012;Stumpf et al., 2013;Bogaard and Greco, 2016).This is why a simple relationship between deformation and rainfall might not be apparent and the detection of this relationship requires displacement and rainfall measurements with high temporal acquisition frequency.
Due to decorrelation of the signal for example in vegetated areas, several advanced processing techniques were developed, the two most common of which are termed Persistent Scatterer interferometry (PS-InSAR, Ferretti et al., 2001;Adam et al., 2004;Hooper et al., 2004) and Small Baseline interferometry (SBAS, Berardino et al., 2002;Schmidt and Bürgmann, 2003;Usai, 2003).
Recent evolutions of post-processing techniques aim to combine information of pixels that contain dominant scatterers, like man-made structures, and pixels that contain mainly natural targets (Hooper, 2008;Ferretti et al., 2011;Hooper et al., 2012;Spaans and Hooper, 2016).
Since PS-InSAR and Small Baseline techniques are optimized for different scattering models, their results can be combined to increase spatial sampling including both groups of pixels with dominant scatterers (i.e.persistent scatterers) and pixels with natural terrain that decorrelate little over short time intervals (Hooper, 2008).Other advanced techniques try to maintain information of both groups of pixels by identifying pixels with similar scattering characteristics (sibling pixels in the following) via a combination of adaptive spatial filters and statistical tests on the amplitudes of SAR scenes in a stack (Ferretti et al., 2011).Because most multitemporal InSAR approaches estimage an average coherence for the whole stack, it is overestimated for some interferograms and underestimated for others.Spaans and Hooper (2016) proposed a technique that identifies sibling pixels in a single master stack by exploiting characteristics of both SAR scenes and interferograms, in order to use only pixels that belong to the same set of siblings during the coherence estimation for individual interferograms.
In the present work, we used DORIS (Kampes and Usai, 1999) and ROI-PAC (Rosen et al., 2004) to process the SAR data from the ASI COSMO SkyMed constellation, while we used GMT-SAR (Sandwell et al., 2011) for the processing of Sentinel 1A.Different works applied successfully standard InSAR approaches with nearly continous spatial coverage, using L-band data in areas that are characterized by the absence of man-made structures (Strozzi et al., 2005;Roering et al., 2009;Handwerger et al., 2013;Tong and Schmidt, 2016).Because of general low coherence in our study area we rely on the Small baseline implementation of the Stanford Method of Persistent scatterers (in the following StaMPS Hooper et al., 2007;Hooper, 2008) to post-process our interferometric results.StaMPS is designed to identify pixels that are coherent in most interferograms, based on the statistical analysis of amplitude and phase characteristics of both the scenes and the interferograms that are used in the survey (see Hooper et al., 2004Hooper et al., , 2007;;Hooper, 2008, for details).
The interferometric phase is not entirely caused by deformation, but contains several terms that can be regarded as noise if deformation measurement is the primary goal of the analysis (e.g.Agram and Simons, 2015).The phase term that is caused by topography can be calculated and subtracted with an external digital elevation model (e.g.Bürgmann et al., 2000).Errors in the elevation model cause a pattern in single interferograms that is proportional to the perpendicular baseline of the interferograms (Massonnet and Feigl, 1998), while in time series approaches a signal is introduced that resembles the baseline history (Fattahi and Amelung, 2013).Although different processing schemes provide the possibility to correct for the phase term due to DEM errors (e.g.Berardino et al., 2002;Hooper et al., 2007;Fattahi and Amelung, 2013), an accurate digital elevation model may aid during pixel selection and can improve the deformation signal (Bayer et al., 2017a).This is why a high resolution digital surface model (2 m) derived from aerial photos acquired in 2008 was used to calculate and correct the effects of topography on the interferometric phase.Another important source of noise is the atmospheric phase delay (Fattahi and Amelung, 2015;Agram and Simons, 2015) and different techniques were developed to correct for the atmospheric phase screen relying either on a correlation between topography and phase or on the calculation of the phase delay from atmospheric models (Bekaert et al., 2015a,b).Because atmospheric disturbances cause patterns in interferograms that are large with respect to a signal that is caused by landsliding, an efficient way of reducing noise in the InSAR time series is the selection of a stable reference area as close as possible to the deforming feature under analysis (Casu et al., 2006).
Interferograms cover large areas from river catchments (e.g.CSK) to whole mountain ranges (Se1), while single landslides produce local deformation signals with respect to the whole scene (e.g.Roering et al., 2009;Colesanti and Wasowski, 2006;Handwerger et al., 2013;Wasowski and Bovenga, 2014).Since the landslide density in our study area is high (on average > 20 %), we are interested in both large scale information, like the detection of landslides, as well as local aspects of the signal, like deformation time series for different landslides.We first performed the interferometric processing on the whole study area in order to identify clear InSAR signals that can be related to landsliding.We considered a line of sight signal to be attributed to a landslide if a consistent spatial difference in the mean velocity is observed that does not cross rivers or ridges (Roering et al., 2009) and if landslides were mapped by the regional inventory.Then, we reprocessed local subsets of the interferograms that include clear landslide signals or landslides that were known to be active, and calibrated the processing parameters to suit the specific landslide under analysis.In particular we were careful in chosing the stable reference area and we adjusted filter dimensions and unwrapping parameters.Because DEM errors or false DEM error estimates may propagate into the temporal displacement signal as a baseline correlated scatter (Fattahi and Amelung, 2013;Bayer et al., 2017a), we paid attention to the correct modelling of the DEMerror term.Two of our InSAR datasets were used in a previous study to assess landslide motion that was induced the excavation of a double road tunnel and the InSAR measurements were in good agreement with traditional inclinometer measurements indicating milimeter accuracy for the periods of high sampling frequency (Bayer et al., 2017b).

Combination of interferometric datasets
COSMO SkyMed acquired SAR imagery on two descending tracks and one ascending orbit (Fig. 2 a and b), while in case of Sentinel we processed one acending and two descending Sentinel tracks over our study area (Fig. 2 a and c).It is possible to project line of sight InSAR measurements along the downslope direction given slope and aspect information from a digital elevation model (e.g.Hilley et al., 2004;Handwerger et al., 2013Handwerger et al., , 2015)).This is not only an intuitive way of representing InSAR measurements for landslide applications, but it delivers also a common coordinate system and where interferograms from different orbits sensed the same area with different viewing angles, we combined interferograms from all available tracks by backprojecting them from the line of sight onto the downslope direction.If more than one InSAR dataset is covering a given area of interest, we choose pixels that are common to all dataset.Similar approaches that involve the resampling from one geometry into the other were used to decompose the signal into vertical and horizontal components along the east-west direction (e.g.Tofani et al., 2013;Eriksen et al., 2017) and, as in our approach, are associated to a loss in detail.In order to solve for deformation in the single time steps with a least squares inversion, a closed interferogram network is required, where each scene of each interferogram is either a master or a slave in another interferogram (e.g.Schmidt and Bürgmann, 2003;Hooper, 2008).To achieve that, we modified one date at the beginning of each dataset, in order to have one common date in all datasets and form a closed interferogram network.The date was chosen in order to minimize the period (less than 10 days) between the modified date and the original dates in the different datasets.Then we solved for the displacements in the single timesteps using the inversion implemented in StaMPS (Hooper, 2008).
This approach yields a series of combined downslope projected InSAR datasets.Their details are listed in Tab. 1 and their spatial coverage is illustrated in Fig. 2 b and c.
One advantage of combined InSAR datasets is that the sampling frequency is higher and that the timesteps are covered by a higher number of observations.For example, during the period of 2013 the COSMO SkyMed datasets have sampling frequencies that range from 16 to 48 days in the different orbits.After the combination of the interferograms the average acquisition rate of the combined dataset is 11 days, with values ranging from 1 to 36 days.For Sentinel, the sampling rate is generally 12 days before interferograms are combined and reach on average 8 days after inverting the back-projected interferograms from different orbital tracks.In addition to the aforementioned loss in detail, the presence of high frequency noise in the inverted time series results and a speckle like pattern in space are a common problem.The speckle pattern in the velocity maps derive in our view from the difference in pixel positioning.We used a local regression filter (lowess, Cleveland, 1979;Cleveland and Devlin, 1988)  on the conditions in the adjacent units, we decided to simplify the nested landslides to a more simple larger perimeter.Occasionally deforming areas are located outside the mapped landslide perimeters, in which case we retraced the perimeters of the InSAR detected landslides, based on the topographic map, field observations, a regional digital elevation model, and aerial photography.
InSAR displacement patterns are visible for 13 landslides in pelitic flysch formations and 12 landslides in chaotic clayshales.We analyse the morphometric characteristics of the remapped, InSAR-detected landslides in the two main lithologies (Fig. 4) and found that differences can be observed in particular with respect to width, length, area and slope angles of the landslides.Landslide length varies between 1100m to 3360m in pelitic flysch formations (mean: 2137m), while landslides in the chaotic clay shales can be longer with landslide lengths ranging from 750m to 3700m (mean: 2137m, Fig. 4 a).On the other hand, the slope failures in pelitic flysch are typically wider (mean: 1566m, var: 600m to 3600m) than their equivalents in the chaotic clay shales Comparing landslide width and length (Fig. 4 e), it can be stated that landslides in the pelitic flysch units can become wider with respect to their length and as mentioned earlier, are often complex landslides with characteristics of both rotational and translational sliding that may transition to classical earthflows.Since they are often composed by numerous geomorphological units or minor nested landslides and the mechanical properties correspond to a soil rather than a rock, they can often be classified as complex earthslides (Cruden and Varnes, 1996).On the other hand the In-SAR detectable landslides in the clay shales have lower width/length ratios, and in most cases they exhibit a bowl-shaped source area, an elongated transport zone and a lobate toe, which is why they can be classified as earthflows (Cruden and Varnes, 1996;Hungr et al., 2001Hungr et al., , 2012)).similar landslide types may occur in both clay shales and pelitic flysch formations.
We detected no clear InSAR signal on landslides that were smaller than approximately 500 m in width or length.In our experience, deformation of small landslides in our study area do not cause clear spatial signals in interferograms, because only few sparse pixels are selected on structures or exposed soil on the landslides.Often these pixels look noisy in the wrapped interferograms and due to large gaps between these sparse pixels, integer phase cycles might not be resolved correctly by unwrapping algorithms and hence the signal would need verification by ground truth data.

Differences in spatial displacement patterns
We analysed the mean velocities of all pixels that were located inside our selected landslides for two different periods.The first period is covered by COSMO SkyMed (2009to 2015) and includes the aforementioned long enduring rainfall events of the years 2013 and 2014.Our Sentinel dataset covers the period between autumn 2014 and summer 2016.
If all pixels inside mapped landslide perimeters are considered, a certain number of pixels reports low velocities beause they are located in parts of the landslide that are not moving.These pixels become more frequent with increasing landslide area.It is also possible that some pixels are noisy including either false high velocities or false low velocities.Both the stable and the noisy pixels bias the density distribution of velocities.To overcome this problem, we decided to use the 75th percentile of the inferred downslope velocity (Fig. 4 f) on a single phenomena (maximum velocity in the following), because it gives less weight to pixels in stable areas and excludes most noisy pixels with false high velocities, being a more robust estimater for landslide velocity.
During the year 2009-2014 (COSMO SkyMed) numerous pixels with mean velocities ranging from 20 to 40 mm/a were registered on the complex landslides of the pelitic flysch formations and only single pixels report mean velocities higher than 50 mm/a (Fig. 4 f).On the landslides in the clay shales, pixels with lower mean velocities prevail and again only single pixels move faster than 50 mm/a.During the Sentinel dataset (2014-2016), the shape of the distribution changes and the differences between the landslide types become less articulated (Fig. 4 f).In the CSK dataset (2009-2014) maximum velocities of pelitic flysch earthslides are higher than those registered on clay-shale earthflows and slowed down in the Sentinel dataset (2014)(2015)(2016).Maximum earthflow velocities remained more stationary throughout all datasets, although velocities and variance are slightly higher compared to the pelitic flysch in the Se1 dataset (Fig. 4 f).Because decorrelation (Zebker and Villasenor, 1992) was a major problem in our X-band dataset and unwrapping was always more challenging compared to Se1, we interpret this difference as a techincal advantage of a C-band radar with high acquisition frequency in a study area with sparse vegetation and few man-made structures and outcrops.

Local variations of different landslide types
Due to regional differences in mechanical characteristics of host rock, landslide material as well as local variations in rainfall, landslides may display a large range of kinematics.Here, we document eight selected cases, four of which involve pelitic flysch formations and four involve the chaotic clay shales.In order to exclude the possibility that our temporal deformation patterns are biased by our InSAR dataset, we chose these cases in different parts of our study area that in the part of the landslide that is known to be active and the signal in the different dataset is measured on the same physical objects on the ground.On the other hand, the cases of Creda and Corzago/Lama di Monchio show that Sentinel selects more pixels in areas with sparse man-made structures and moderate vegetation and can better resolve locally high displacement rates.Because the CSK dataset has a shorter wavelength and a lower acquistion frequency compared to Se1, it is more affected by temporal decorrelation (Zebker and Villasenor, 1992).This probably why for CSK less pixels are selected in the strong deforming area and the higher velocities in Se1 are not caused by an acceleration of the landslide bodies, but by the technical differences between CSK and Sentinel.In our view it is possible that here displacement rates were even higher during 2013 (CSK), but were not fully resolved by the X-band satellite.
The four representative earthflows are reported in Fig. 6.Vidiciatico (Fig. 6 a) and Gaggio Mon-tano (Fig. 6 b) are medium sized phenomena, while Monte Ombraro (Fig. 6 c) and Sestola (Fig. 6 d) are large earthflows.All of them exhibit the distinctive morphological features of earthflows with characteristic source area and accumulation zones (Simoni et al., 2013).For the sake of this work (Fig. 6), we mapped the entire earthflow complex as a whole although multiple coalescent earthflows can be further recognized and separately mapped on a geomorphological basis.At Vidiciatico and Gaggio Montano, the signal appears consistent throughout all datasets.In the case of Gaggio Montano, slow displacement rates were measurable also in older InSAR datasets (Corsini et al., 2006;Villi et al., 2016) and confirm the sustained slow movement typical of these phenomena.
The comparison between COSMO SkyMed and Sentinel datasets indicate that locally higher mean velocities are measured with the C-band satellite.This difference is well visible in the central part of Monte Ombraro (rectangles in Fig. 6 c).

Deformation responses to long enduring rainfall events
Landslides are complex hydrological systems and, as a consequence, a simple relationship between rainfall and displacements is rarely evident (e.g.Cervi et al., 2012;Berti et al., 2012;Bogaard and Greco, 2016).Still, rainfall, infiltration and the resulting variation in pore pressures can be considered the main drivers of landslide motion (e.g.Iverson, 2000;Berti and Simoni, 2012;Berti et al., 2012;Handwerger et al., 2013).In the following, we document the velocity time series of the highlighted pixels in Fig. 5 and Fig. 6 together with rainfall derived from rain gauges that are distributed over our study area (Fig. 3 a).
We report the velocity time series during 2013, because the sampling frequency was high for all CSK datasets (Fig. 7 a and b) and 2015 (Fig. 7 d and e), because all cases were covered completely by all Sentinel datasets.The rainfall data from all rain gauges was averaged for each day in order to obtain a mean value for the study area.Since our combined displacement measurements have average temporal sampling rates between 8 (Se1A) and 12 (CSK) days, as well as maximum time spans between acquisition of up to 36 days, we compared 20 days cumulative rainfall to our displacement measurements (black lines in Fig. 7 c and f).We also compared displacement rates to rainfall events that following Berti et al. (2012) were defined as a period of continous rainfall separated from the next one by period of at least 3 days with rainfall below 5 mm (grey bars in During the year 2013 persistent rainfalls started by the end of January (first black bar in Fig. 7 c) which increased displacement rates only at La Borra and Camugnano, while the velocities on the other landslides did not increase (Fig. 7 a).March and April 2013 were characterized by sustained rainfalls spanning 35 days, which define a single event characterized by moderate intensity (ca.8 mm/day) but high absolute rainfall (ca.270 mm, second black bar in Fig. 7 c).This rainfall event caused notable accelerations on most landslides of the pelitic flysch, and is deemed the triggering rainfall for the landslide reported by Berti et al. (2017).The first landslide to respond to these rainfalls was La Borra, followed by Camugnano and Creda, while Corzago/Lama di Monchio did not show any increase in displacement rate.It is worth to recall that the discussed landslides are spread over a relatively large territory (3600 km 2 ) and that they are covered by different datasets.Still, the observed InSAR signals are coherent and often synchronized, which illustrates well the relationship between acceleration and meteorological forcing.For all landslides that show an increase in displacement rate, peak velocities are reached between 20 and 60 days after the onset of the long rainfall event in March/April.The lag to peak velocity is consistent with the observations made by Handwerger et al. (2013) on slow moving earthflows in mélange type rocks of Northern California.We believe that this lag is not surprising for deep seated landslides, because different processes like surface infiltration and pressure diffusion (Iverson, 2000), as well as the saturation of fissures (Krzeminska et al., 2012) can delay the increase in pore pressure.
During the spring 2015, total recorded rainfall was lower than the two previous years and also single rainfall events never reached magnitudes similar to the April 2013 event.However, the trend of 20-days cumulative rainfall shows that the wet season started in mid January with major rainfall events in January and February (Fig. 7 f).Cumulated rainfall remained high until May, although rainfall was interrupted by several dry days.The displacement signals registered on   2012), while the black line depicts daily rainfall that takes into account the preceding 20 days of rainfall.We chose 20 days because this timespan is in good agreement with our definition of rainfall events.d) Displacement rates derived from Sentinel for the year 2015 on rockslides in pelitic flysch.e) Velocity time series for earthflows in mélange type rocks.f) Daily rainfall that takes into account the preceding 20 days and it is plotted together with rainfall events for the year 2015.
landslides confirm the lithology-dependent variability described above.Higher peak velocities and accelerations are found on pelitic flysch landslides while the seasonal displacement rates resemble smoothly the pattern of the cumulated 20 days rainfall (black line in Fig. 7 f) for all landslides in clay shales.Those accelerating during spring 2015, reacted also after the rainfall events in January and February.If we consider the January event as the onset of the spring rainfalls, the lag to peak landslide velocity is 20 to 60 days also in early 2015.

Relationship between deformation response and rainfall
Mechanical aspects, like the tendency of the sheared clay material to contract or dilate (Iverson, 2005) or hydrological aspects like priamary diffusivity (Iverson, 2000), the development of effective drainage pathways in time (Bogaard and Greco, 2016) impose non linearities on the relationships between rainfall, pore pressure and displacement.Different types of fissures may develop in a landslide deposit (Stumpf et al., 2013) and may play a major role for the build up of porepressure (Krzeminska et al., 2012;Bogaard and Greco, 2016).Corominas et al. (2005) showed that the relationship between displacement rate and pore pressure is not linear and can be approximated by a third degree polynomial.Massey et al. (2013) have shown that a high correlation between rainfall and pore pressure exists only if more than 30 days of preceding rainfall is taken into account.
Hence, whenever precipitation is compared to deformation or landslide failure in general, a period of antecedent rainfall or a definition of rainfall event is necessary, although values may vary locally (see for example Guzzetti et al., 2007, for a compilation of regional rainfall thresholds).Tong and Schmidt (2016) binned the precipitation that occurred in each timestep of InSAR measurements to compare deformation with precipitation, while Corsini and Mulas (2016) reported that, adding the antecedent 40 days of rainfall to each rainfall cause high correlation between displacement measurements from total stations and rainfall on a Northern Apennines rock slide.
Chosing the right period of antecedent rainfall is often subjective, which is why we tested the relationship between acceleration/velocity and rainfalls for different periods of antecedent rainfall (Fig. 8).Also the representation of deformation affects the degree and signicance of the correlation.We chose to compare velocities/displacement rates as well as a change in displacement rate (acceleration), which can be defined as the actual response to rainfall.In analogy to Handwerger et al. (2013), we refer to the period between the beginning of a rainfall event and the timestep during which an acceleration occured on a landslide as response time.
A comparison between displacement rate and rainfall on landslides in pelitic flysch formations yielded only moderate but statistically significant correlation coefficients (> 0.25) with p-values well below 0.05, if more than 30 days of preceding rainfall is summed to the daily values and remains high if longer periods of preceding rainfall are considered (Fig. 8 a).Massey et al. (2013) reported similar high correlation between pore pressures and rainfall for long periods of preceding rainfall on a complex landslide in marine sand and siltstones in New Zealand.Pore pressure decay can take longer than its rise (Berti and Simoni, 2010) and hence periods of relatively high pore pressures can extend into the dry season, especially if motion started late during the wet season.
In analogy, once motion is triggered, periods of elevated displacement rates can be long and may extend into the beginning of the dry season.
Moderate and statistically signficant correlation coefficients (> 0.25) with P-values well below 0.05 were obtained for the relationship between rainfall and acceleration if more than 20 and less than 60 days of preceding rainfall were considered (Fig. 8 b).This is consistent with InSAR derived seasonal response times of slow moving earthflows in Northern California (Handwerger et al., 2013) and a rockslide in the Northern Apennines (Corsini and Mulas, 2016).On the contrary, periods of acceleration and deceleration are short and if compared to a long window of antecedent rainfall, periods with no acceleration or deceleration will often have high total rainfalls and as a consequence low correlation-coefficients.
A comparison between displacement rate and rainfall for landslides in clay shales yielded lower correlations (< 0.25) with P-values oscillating around 0.05, if more than 30 days and less than 80 days of preceding rainfall were considered, suggesting a longer response time (Fig. 8 c).No clear relationship between acceleration and deceleration can be observed for any period of anteceding rainfall (Fig. 8 d).We observe that, in analogy to periods of acceleration for flysch landslides, periods of high displacement rates were short and corresponded to the spring rainfall 2013 and 2015.On the other hand, acceleration and deceleration were longer and low in magnitude, which is why no correlation exists.
Since peak correlation coefficients with low p-values were measured for periods between 20 and Velocity (mm/month) Precipitation (mm) Precipitation (mm)   60 days of preceding rainfall, we compare in the following the velocities and accelerations in the two lithotypes to daily rainfall that takes into account the preceding 40 days of rainfall (Fig. 9).
Landslides in pelitic flysch show a weak relationship between velocity and rainfall (Fig. 9 a).
For landslides in chaotic clay shales a moderate relationship between velocity and rainfall can be observed (Fig. 9 c), while no correlation was measurable between rainfall and acceleration.Still, also for the landslides in clay shales periods of high velocities and acceleration were measured during the wet season.

Discussion
Slow moving deep-seated landslides and earthflows can display a wide spectrum of displacement styles that include steady linear displacement, moderate seasonal variation in deformation rates, but they may also transition to catastrophic failure (Keefer and Johnson, 1983;Petley and Allison, 1997;Petley et al., 2005).Only few studies exist that investigated temporal displacement patterns on the catchment scale (Roering et al., 2009;Handwerger et al., 2013Handwerger et al., , 2015;;Bennett et al., 2016;Stumpf et al., 2017).In the present study the InSAR derived displacement rates remain typically below 10 cm/year in the downslope direction and only large landslides delivered clear signals.
Detailed seasonal variation and linear displacement rates were documented for two years, although deformation maps derived from COSMO SkyMed and Sentinel datasets spanning more than five years suggest that these slow displacements are persistent in time.We would like to point out that these deformations are consistent with seasonal variations in long enduring creeping phases (Petley and Allison, 1997).Several reactivations that were registered during the spring rainfalls of 2013 were too fast to be measured by InSAR.Berti et al. (2017) described in detail a 30 m deep failure that occurred on 6th of April in pelitic flysch rocks.Total displacements of up to 45 meters were inferred from differencing pre-and post-failure topography and eye witnesses report displacement rates of ca. 10 m/h.
The CSK dataset (2009-2015) allowed us to measure higher velocities on pelitic flysch landslides, compared to the earthflows in the clay shales in particular during the year 2013.Independent from landslide type and bedrock lithology, overall displacement rates are lower in the Sentinel dataset (2015-16) compared to COSMO-SkyMed.However, locally higher displacements rates were resolved by Sentinel due to higher acquisition frequency and good coherence in areas with sparse man made structures.We combined interferograms from different viewing geometries by choosing common pixels and backprojecting the data along the downslope direction before solving for displacement time series.Although this approach reduced the detail of our datasets, it helped to better constraint periods of acceleration and deceleration because more interferograms span the same epic.
The higher mean velocities during the CSK period are in line with reports, that state a high number of reactivations of landslides in pelitic flysch units and only few reactivations of earthflows in clay shales (Pizziolo et al., 2013(Pizziolo et al., , 2015)).Because only a reduced number of the mapped landslides in the study area showed high displacements during our observation period and an even smaller number of landslides were measurable by means of radar interferometry, we analysed only a reduced sample, compared to the the total area that is classified as landslide affected terrain.However, the InSAR derived displacement time series confirm that in 2013 landslides in pelitic flysch rocks accelerated faster and reached higher peak velocities, while displacement responses remained lower in magnitude on the earthflows in clay shales.The fact that the recorded rainfall events during the observation period consisted of sustained rainfall reaching only moderate intensities, suggests that these rainfalls had a critical impact on the landslides in pelitic flysch.Handwerger et al. (2015) proposed that shallower landslides are more sensitive to variations in rainfall, as porepressure at the slip surface rises faster and Bennett et al. (2016) showed that deeper landslides responded later to long term drought conditions.Different landslide depth might in part explain the different response times in the same lithological formations, but it is unlikely to be the root of the different kinematics of the landslides in the pelitic flysch compared to the ones in chaotic clay shales.Instead it can be hypothesized that the specific rainfall pattern during spring 2013 had different influence on the pore-water pressure regime of the two lithotypes, because they have in all likelihood signifcantly different hydrological characteristics.Pelitic flysch units have higher hydraulic conductivity due to fissures and discontinuities that tend to close with depth but also to open with deformation.The water table can be located at various depth from the surface and exhibit significant seasonal variations (Cervi et al., 2012;Borgatti et al., 2006).In our view large volumes of infiltrating water are required to raise the pore-water pressure within the pelitic flysch to trigger deformation.On the other hand, clay-shales have lower hydraulic conductivities and although networks of fissures may also be present (Berti and Simoni, 2010), they may not persist at depth where pore-water pressures remain high (close to hydrostatic) throughout the year Berti et al. (2013).It is possible that failures in clay-shales require rainfall with high intensities rather than long durations.
The direct statistical comparison between precipitation and deformation is complicated by additional factors.On the one hand it is sensitive to the period of antecedent rainfall that is taken into account.On the other hand, the representation of deformation makes a difference, because in pelitic flysch formations acceleration and hence the period during which displacement rate changes appears to be stronger correlated with precipitation compared to velocities.A possible explanation for the low correlation between precipitation and velocity is that, once displacement is triggered by the increase of pore-water pressure, velocity is regulated by the friction characteristics of the material, landslide geometry and inertia, together with the decay of pore-water pressure that is slower than its rise (Berti and Simoni, 2010).Acceleration appears more clearly related to rainfall (Fig. 9 b) and in all cases, periods of high acceleration were measured during the wet seasons in autumn and spring, while periods of deceleration occurred mostly during the dry summer.Both visual analysis of the displacement time series and statistical comparison between acceleration and rainfall suggests that landslide response times for our landslides vary between 20 and 60 days, which is consistent with the findings of Handwerger et al. (2013).Corominas et al. (2005) demonstrated that the relation between pore pressure and displacement is complex for periods shorter than 20 days and displacement rates higher than 60 mm/month, which were however not resolvable by our InSAR dataset.Here, we observed a weak but linear relationship between rainfall and slow velocities if more than 20 days of preceding rainfall were taken into account.A possible explanation for the linear nature of this relationship is the fact that only slow deformations were registered and that the aforementioned non-linearities become more important at higher displacement rates and sudden failures, which can not be measured by radar interferometry.

Conclusions
In the present work we document InSAR derived spatio-temporal displacement patterns of 25 landslides in three river catchments of the Northern Apennines, where pelitic flysch rocks and chaotic clay shale formations constitute most of the substrate.We show that different lithologies host different typologies of landslides.The complex landslides of the pelitic flysch units responded with abrupt acceleration and higher displacement rates to long enduring rainfalls.Chaotic clay shales host classical earthflows that responded with overall lower displacement rates and smoother acceleratons.The difference in kinematic behaviour can be explained by different hydrological and geomorphological characteristics.In particular permeability and slope angle, of bedrock and landslide material might play an important role.The landslides where creeping rate was related to rainfall, responded between 20 to 60 days after the onset of the spring rainfall.Because displacement rates were generally low, the relationship between different displacement measures and rainfall is approximately linear if longer periods of preceding rainfall was taken into account.

Figure 2
Figure 2: a) The area of interest is located in the Northern Apennines of Italy and is covered by several SAR frames in both ascending and descending viewing geometries.b) For COSMO SkyMed one ascending and two descending swaths were available.A small area is covered by all three datasets (CSK3), while most parts of our study area are covered by either two (CSK2a and CSK2b) or one dataset (CSK1).c) Also for Sentinel 1A one ascending and two descending tracks were processed, which results in a larger area covered by either two (S2a and S2b) or three InSAR datasets (S3).Details regarding the different datasets and their combination can be found in Tab. 1.

Figure 3
Figure 3: a) Calibrated InSAR results for all landslides that caused a clear deformation signal in the COSMO SkyMed datasets.Black stars indicate rain gages that were used to retrieve regional rainfall that will be discussed together with displacement signals in section 4.2.b) Calibrated InSAR results for Sentinel 1A.The labels in a) and b) are localities, where distinct landslides were detected.From West to East and from South to North we have: Bs: Boschi di Valoria, Fr: Frassinoro, Cs: Casolare, Ro: Roncolo, Mf: Montefiorino, Lm: Lama di Monchio, Lb: La Borra, Se: Sestola, Vi: Vidiciatico, Sg: San Giacomo, Gm: Gaggio Montano, Ct: Castelluccio, Bc: Borgo Cappanne Mo: Monte Ombraro, Sv: Suviana, Ca: Camugnano, Cp: Castiglione dei Pepoli, Cr: Creda, Bg: Il Borgo.Red highlighted labels indicate selected landslide cases that are discussed in detail in section 5. Blue labels are cases that were only visible in the Sentinel dataset.

(
mean: 567m, var: 196m to 1450m Fig.4 b).Due to their larger variance in width, pelitic flysch landslides cover larger areas (mean: 3.06km 2 , var: 0.16km 2 to 10km 2 ) compared to the those in clay shales (mean: 1.2km 2 , var: 0.12km 2 to 4.34km 2 , Fig.4 c).Slope angles on landslide prone terrain in chaotic clay shales reach on average 11 degrees (var: 8 to 15.5), while slopes in pelitic flysch formation maintain on average angles around 12.4degrees (var: 9.5 to 14.9 , Fig.4 d, Fig.1c in supplementary I).If we use slope angle as a proxy for mechanical properties of both bedrock and regolith, we can infer that chaotic clay-shales have lower shear resistance than pelitic flysch.

Fig. 4 eFigure 4 :
Figure 4: Probability density distributions calculated with a gaussian kernel of a) landslide lengths, b) widths, c) area, d) slope angles and e) the relationship between width and length, normalized by the landslide area for InSAR detectable landlsides in our study area for the two main lithologies.f) 75 percentile of the mean downslope velocities measured on the landlslides in pelitic flysch formations and chaotic clay shales.

Figure 5 :Figure 6 :
Figure 5: Deformation maps for selected rockslides in pelitic flysch formations.CSK and Se1 results a) at La Borra, b) at Creda, c) at Corzago/Lama di Monchio, d) at Camugnano.Due to the high spatial resolution, the number of registered pixels are higher for COSMO SkyMed, yet Sentinel selects also points in areas, where COSMO SkyMed was decorrelated.Displacement patterns vary in space for the different periods and the highlighted pixels were used to produce the time series in Fig. 7 a) and d)

Fig. 7 c
Fig. 7 c and f).First, peak velocities are higher for most pelitic flysch landslides during Spring 2013, reaching values up to 30 mm per month in the downslope direction (Fig. 7 a), while peak

Figure 7
Figure 7: a) Velocity time series for landslides in pelitic flysch during for the year 2013.b) Velocity time series for earthflows in the chaotic clay shales for 2013.c) The grey bars are rainfall events that were defined by the methodology proposed by Berti et al. (2012), while the black line depicts daily rainfall that takes into account the

Figure 8 :
Figure 8: Correlation coefficients and P-values for the relationships between a) velocity and rainfall and b) acceleration on landslides in pelitic flysch formations and different days of preceding rainfall added to the daily rainfall values.c) Correlation coefficients and P-values for the relationship between velocity and d) acceleration on landslides in chaotic clay shales.

Figure 9
Figure 9: a) Relationship between velocity on landslides in pelitic flysch and rainfall if 40 days of preceding rainfall are summed to the daily values.b) Relationship between acceleration and rainfall on pelitic flysch landslides.c) Comparison between velocity and d) acceleration of landslides in chaotic clay shales and rainfall.

Table 1 :
Detailed information about single and combined InSAR datasets.Ifgs. is an abbreviation for interferograms