Establishing the timings of individual rainfall-triggered landslides using Sentinel-1 satellite radar data

Heavy rainfall events in mountainous areas can trigger thousands of destructive landslides, which pose a risk to people and infrastructure and significantly affect the landscape. Landslide locations are typically mapped using optical satellite imagery, but in some regions their timings are often poorly constrained due to persistent cloud cover. Physical and empirical models that provide insights on the processes behind the triggered landsliding require information on both the spatial extent and timing of landslides. Here we demonstrate that Sentinel-1 SAR amplitude time series can be used to constrain landslide timing 5 to within a few days and present three methods to accomplish this based on time series of: (i) the difference in amplitude between the landslide and its surroundings, (ii) the spatial variability of amplitude between pixels within the landslide, and (iii) geometric shadows cast within the landslide. We test these methods on three inventories of landslides of known timing, covering various settings and triggers, and demonstrate that, when used in combination, our methods allow 20% of landslides to be timed with an accuracy of 80%. This will allow multi-temporal landslide inventories to be generated for long rainfall 10 events such as the Indian summer monsoon, which triggers large numbers of landslides every year and has until now been limited to annual-scale analysis.

In order to obtain timed landslide information for inventories associated with long or successive rainfall events, we propose a two-step process, whereby landslide locations are mapped as polygons using optical satellite imagery, and the timings of individual landslides are then obtained from SAR time series. In this paper we address the second of these steps. We use 60 Sentinel-1 time series over inventories of landslides whose timings are already known to test three potential landslide timing methods individually and in combination.

Case study events
We used three published polygon inventories of landslides whose timings are known a-priori to test and develop landslide timing methods. We filtered each inventory to remove landslides smaller than 2000 m 2 , so that each landslide was expected to timings before, during and after a defined "co-event" window of 6 months relative to the real event timing.
valleys we consider here, only Bhote Kosi was close enough to the epicentre to be affected by that event (Martha et al., 2017). 90 Since the co-event pair of SAR images for Bhote Kosi (24 April -18 May 2015) spans both the Gorkha earthquake on 25 April and the Dolakha aftershock on 12 May, these two earthquakes can be considered as a single triggering event in Bhote Kosi.

Theory: SAR backscatter and landslides
A SAR satellite actively illuminates the Earth's surface with microwave energy, and records the phase and amplitude of the returned signal. The difference in phase between two images acquired over the same area at different times can be used to track 95 the movement of the Earth's surface, for example movement on a fault during an earthquake, while the amplitude describes the strength of the backscattered SAR signal. The power of the signal transmitted P r and received P t by the sensor are described by Eq. 1, where λ is the wavelength, G 2 is the two-way antenna gain and R is the slant range (Small et al., 2004).
This equation is solved to obtain x 0 , the backscatter coefficient, which can be either σ 0 , γ 0 or β 0 depending on whether the 100 integration is carried out in the ground (ellipsoid) plane, the plane perpendicular to the look direction or the slant-range plane respectively (Small et al., 2004). Different studies have demonstrated that all three of these backscatter coefficients can be applied to detect vegetation removal due to landslides and other processes such as deforestation and wildfires (e.g. Ban et al., 2020;Belenguer-Plomer et al., 2019;Bouvet et al., 2018;Esposito et al., 2020;Hernandez et al., 2021;Konishi and Suga, 2018;Mondini, 2017;Mondini et al., 2019;Motohka et al., 2014). Here we used γ 0 .
105 SAR backscatter is dependent on a number of factors, including the polarisation and wavelength used by the SAR system, the local slope orientation relative to the SAR sensor and the roughness and dielectric properties (e.g. soil moisture, presence of vegetation) of the material that the microwave energy interacts with at the Earth's surface. Sentinel-1 acquires C-band SAR data with a wavelength around 5.5 cm in two polarisations: "VV" (vertical polarisation) and "VH" (cross polarisation). We tested both of these polarisations, but found VV to perform better than VH so present only the results for VV. VV data have 110 also been acquired more consistently throughout the lifetime of Sentinel-1 than VH. In general, for vertically polarised SAR images, rougher surfaces result in increased backscatter, as does increased soil moisture. However, the relationship between these properties and the SAR amplitude is not simple: roughness has a stronger effect in locations with a high incidence angle (Baghdadi et al., 2016;Dubois et al., 1995), while changes in soil moisture have a larger effect at low incidence angles (Baghdadi et al., 2016).

115
Landslides alter the local topography (and therefore the local incidence angle) of the landscape through the movement of material and remove vegetation, which alters the dielectric properties and roughness of the Earth's surface. For this reason, landslides can result in both increases and decreases in amplitude. In fact within a single landslide, the amplitude of some pixels may increase while some decrease (e.g. Tozang landslide, Mondini et al., 2021, Fig 4).

120
To construct our SAR amplitude time series, we used the Google Earth Engine Sentinel-1 ground range detected (GRD) data set. These data are preprocessed following the workflow of Filipponi (2019) to obtain the σ 0 backscatter coefficient in geographic coordinates at a resolution of 20 x 22 m and a pixel size of 10 x 10 m. We then applied the module of Vollrath et al.
(2020) using the SRTM 30 m digital elevation model (DEM) to carry out an angular radiometric slope correction based on the volume scattering model of Hoekman and Reiche (2015). This has the effect of converting from σ 0 (normalised in the ellipsoid 125 plane) to γ 0 (normalised in the plane perpendicular to local satellite look direction). The aim of this step is to reduce the effects of topography on the SAR backscatter. In preliminary testing, we found that γ 0 performed better than σ 0 . The module of Vollrath et al. (2020) also provides a shadow and layover mask that can be used to remove areas that are not well-imaged by the satellite due to the viewing angle and local topography. This masking step is important for landslide studies as they are likely to be carried out in areas of steep topography. For each of our three events, we defined a "co-event" period of approximately six months. We also defined a three-month pre-event period and two-month post-event period immediately before and after the co-event window. These pre-event and post-event image stacks are required in some of the methods outlined in Sect. 2.
The dates that make up the pre-event, "co-event" and post-event time series for each case study are shown in Fig. 1d.
The length of the "co-event" period was defined as 6 months for the Hiroshima and Zimbabwe events. For the three Nepal 135 inventories, this was reduced to 5 months in order to allow a sufficient pre-event images to be acquired following the satellite launch in 2014 and sufficient post-event images to be acquired before the end of July, since few Sentinel-1 images are available can also result in amplitude change. In particular, the rainfall that triggers the landslides will alter the soil moisture content and so may also alter the amplitude of the returned signal. To overcome this, we calculate a background amplitude signal for each landslide. First, we calculated a buffer region between 30 and 500 m around each landslide (Fig. 2a). Then we filtered this buffer to remove any pixels that lie within other landslide polygons and pixels that are dissimilar to those within the landslide, for example pixels located on the opposite side of a ridge, in a river or with different surface cover. In order 150 to assess pixel similarity we calculated three surfaces from pre-event satellite imagery. First, we calculated the normalised difference vegetation index (NDVI) from a single pre-event Sentinel-2 (or, where this was unavailable, Landsat-8) image for each event. Pixels with similar NDVI values are expected to have similar land-cover. Second, we used a stack of N pre-event SAR images i ( Fig. 1) to calculate the mean amplitude A mean,j (Eq. 2) and amplitude variability ∆A mean,j (Eq. 3) for every pixel j through time. Pre-event amplitude and amplitude variability have previously been used by Spaans and Hooper (2016) 155 to identify statistically similar pixels in SAR images.
∆A mean,j = 1 (c,d) Example time series for a single landslide from the Hiroshima dataset using SAR data from Sentinel-1 track 090D. Blue bar shows the duration of the rainfall event during which the landslide was triggered. (c) The median SAR amplitude for the landslide, the local background signal and areas of geometric shadow. (d) SAR amplitude standard deviation for pixels within the landslide For every landslide, we calculated the range of NDVI, A mean,j and ∆A mean,j values seen within the polygon. The buffer of this landslide was then filtered to remove pixels with values outside these ranges. From this filtered buffer, we then calculated 160 the median background amplitude for every image in the co-event time series. A step change in the difference between the median landslide amplitude and the median background amplitude is then used as an indicator of landslide timing. As previously described, landslides can result in both increases and decreases in SAR amplitude. When combining methods, we found that using the increase and decrease in amplitude as separate inputs resulted in better performance than combining these into a single input, for example based on the absolute change in amplitude. Therefore, for this method, we accept both a step increase 165 and a step decrease as a potential indicator of landslide timing.  Ban et al. (2020) observed that in forested and grassland areas, the removal of vegetation due to forest fires led to an increase in the variability of vertically polarised Sentinel-1 γ 0 between neighbouring pixels. Since landslides result in a similar denudation of vegetated areas, we expect that similar effects may occur. Therefore, we calculated the standard deviation of γ 0 within each 170 landslide polygon and used a step increase in this as a potential indicator of landslide timing (e.g. Fig.2d).

Method 3: Geometric shadows
Since SAR is acquired obliquely (with an ellipsoid incidence angle of 31-44 • for the data used here), steep changes in scatterer surface height can result in geometric shadows. The wavelength of Sentinel-1 means that is primarily scattered from the canopy in forested areas, which means that shadows can be cast at the edges of deforested areas if these edges run approximately 175 perpendicular to the satellite look direction (Fig. 2b). Bouvet et al. (2018) developed a method for automatically detecting deforested areas based on these geometric shadows. Since landslides remove vegetation, we expect that shadows should also be cast at the edges of landslides, and that the appearance of new shadows could be used as an indicator of landslide timing.
Furthermore, the three-dimensional shape of the landslide could result in shadows cast within the landslide itself, for example if the landslide has a steep scar. This effect has previously been observed within a large landslide in Nepal by Ao et al. (2020).

180
It is worth noting that, while Bouvet et al. (2018) applied their methods in areas of gentle slopes, the area of a shadow cast by an object of a given height is dependent on slope and aspect: trees of the same height will cast a larger shadow on slopes facing away from the sensor than on those facing towards it. Therefore, we expect this method to be more successful for slopes that face away from the sensor.
In comparisons of multiple inventories of the same event prepared by different people or groups, there are often small 185 discrepancies in the exact size, shape and location of each landslide (Milledge et al., 2021;Pokharel et al., 2021). Since shadow pixels are most likely to lie at the edges of the landslide polygons, it is important not to exclude the edge of a landslide from the analysis. Therefore we extended the area covered by each landslide polygon by 10 m (one SAR pixel) where this did not lead to intersection with another landslide in the inventory. We then identified pixels whose amplitude decreased within this enlarged polygon as shadows. Bouvet et al. (2018) identified shadow pixels as those whose γ 0 value decreased by >= 4.5 190 dB during the deforestation event. We tested values between 3 and 6 dB and also found that a threshold of 4.5 dB performed best. We calculated the mean γ 0 value for every pixel from the pre-event and post-event image stacks and assigned those that decreased by >= 4.5 dB as shadow pixels. The co-event time series of these shadow pixels was then analysed and a step decrease in the median shadow γ 0 relative to the median background γ 0 (Sect. 2.4.1) was used as an indicator of landslide timing.  The size of the peak or trough depends on the magnitude of the increase or decrease, the level of noise elsewhere in the time series and the length of the co-event time series (n dates ). A bigger peak or trough for a time series of the same length indicates a larger step change and less noise and is therefore a more reliable indicator of landslide timing. We therefore apply a peak size threshold to remove unreliable landslide timing estimates. To select this threshold for each method, we use the F1-measure, 205 a statistic that combines both precision and recall. This F1-measure was calculated for a range of peak thresholds using the confusion matrix defined in Table 1 (Fig. 3). Based on this, we require a peak of 0.4 × n dates for the Landslide-Background Method, 0.2 × n dates for the Pixel Variability Method and 0.75 ×n dates for the Geometric Shadows Method. We also assessed whether the level of noise in the time series for each metric (estimated from the variability in pre-event and post-event time series), could be used to indicate whether a timing estimate was likely to be correct, but found this to be less reliable than the 210 convolution peak size.

Results
We used each method in Sect. 2 to assign landslide dates for the five case study areas described in Sect. 2.1. Not all landslides are assigned a date by every method, for example if no geometric shadows are cast within the landslide polygons. Table 2 shows the number of landslides assigned a date by each method and the percentage of these dates that were correct in each 215 case. A baseline was calculated from 1 n dates : the percentage of landslides we would expect to be assigned the correct date by chance for a method with no skill. All three methods consistently perform better than this baseline.

Combining methods
To assess the methods in combination, we take whichever date is predicted the most often for every landslide. Since it is not possible for both a step increase and a step decrease in the Landslide-Background Method to predict the same date, the 220 maximum number of times the same date can be predicted is 3. The number of landslides assigned a date by at least 2 methods and by all 3 methods and the percentages of these that are correct is shown in Table 2. The strong reduction in number of timed landslides when going from an individual method to 2 and then 3 methods in combination underlines the fact that the nature Pixel Variability Geometric Shadows of the change in amplitude varies widely between landslides. However, landslides dated by 2 or 3 methods are correctly dated much more often. Across all 8 tracks, 56 landslides are predicted a date by 3 methods of which 52 (93%) are correct. Fig. 4 225 shows the number of times each date pair is predicted by ⩾ 2 methods and by 3 methods.

Combining tracks
Sentinel-1 acquires data on an ascending track, (moving northwards and looking east) and a descending track (moving southwards and looking west). Therefore, we can generate two sets of dates for each landslide inventory (excluding Bhote Kosi and Buri Gandaki where only descending track SAR data were available). This has several advantages. First, landslides that were 230 not assigned a date using data from one track may be better dated by the second, increasing the number of landslides in the inventory for which a date can be assigned. In particular, landslides that are masked due to foreshortening or layover in one track may be better imaged in the other track. This was advantageous in Trishuli, where both ascending and descending tracks were available and steep slopes meant that large numbers of landslides were masked in each case. 474 and 485 landslides were imaged by the ascending and descending tracks respectively out of a possible 650 landslides in Trishuli. However, 516 235 landslides were imaged by at least one of the two tracks, thus coverage was improved by using both. Second, the acquisition dates of the two tracks are slightly offset, so a landslide that is assigned a date by both tracks has a more precise timing. For Table 2. For each case study, the total number of landslides, the number that are masked due to foreshortening or layover in the SAR images and amplitude timing results. For each method and combination of methods, we give the number of landslides assigned a date followed in brackets by the percentage of these assigned dates that are correct (the specificity of each method). Where timings were combined from multiple methods (M) or tracks (T), the number of these is specified in brackets March 2019, improving the precision from 12 days to 7. This more precise date is also more likely to be correct since it is 240 derived from two sets of independent observations of the landslide.
For each event, we used the ascending and descending tracks to generate dates predicted by at least 2 methods on the same track. Out of all the non-masked landslides in each inventory, 23% were assigned a date in Hiroshima, 21% in Zimbabwe and 14% in Trishuli and of these, 80% of the estimated dates in Hiroshima were correct, 73% in Zimbabwe and 81% in Trishuli (Table 2). As well as being dated by 2 methods on the same track, a small number of landslides were dated either by the 3rd 245 method on that track or by at least one method on the other track. These landslides whose dates were assigned based on 3 or more methods ("⩾ 3M" in Table 2) were assigned the correct date more often than those assigned a date only based on 2 methods on one track ("2M, 1T" in Table 2). We also tested the case where landslides were dated based on the same date being selected by one method from each track, but found that this yielded too many incorrect dates to be useful.

250
We assessed the performance of our dating methods as a function of the landslide characteristics, in terms of pre-event vegetation, landslide area, and slope aspect. This may inform future application of our methods.

Vegetation
In order to assess the effect that vegetation cover has on the methods we propose here, we compared the number of correctly timed, incorrectly timed and untimed landslides with different values of pre-event NDVI (Fig. 5 a-c). We took the maximum 255 NDVI value for each pixel in the year preceding the event and used Sentinel-2 data for Zimbabwe and Hiroshima and Landsat 8 for Trishuli. In all three inventories, the majority of mapped landslides occurred in vegetated areas (0.6 < N DV I < 0.8). In all three cases, a landslide in a more vegetated area was more likely to be assigned a date and this date was more likely to be correct.

260
Another factor that could potentially effect the applicability of the methods we present here is landslide area. Fig. 5d-f shows the distribution of landslides against landslide area. In Zimbabwe and Hiroshima, a higher proportion of larger landslides were assigned a date and a higher proportion of these assigned dates were correct, but these effects were not observed in Trishuli.
We limited our testing to landslides whose area was greater than 2000 m 2 . Since our methods rely on landslides containing multiple SAR pixels in order to calculated the statistics such as the standard deviation, there is likely to be a lower limit on the 265 area of landslides that can be timed that was not reached here.

Aspect
The effect of aspect on landslide timing ability is more complicated than that of vegetation and area, since it is likely to vary between the ascending and descending track SAR. Therefore, in Fig. 6 how likely a landslide is to be assigned the correct time. For the Geometric Shadows Method, a higher proportion of landslides are assigned a date (and therefore exhibit a shadow) on slopes facing away from the sensor. This was expected since the same 275 height difference will cast a larger shadow on a slope facing away from the sensor than one facing towards it (Bouvet et al., 2018). Dates assigned by the Geometric Shadows Method also appear more likely to be correct for slopes facing away from the sensor on Z072A, but this pattern is less clear on Z079D.

Possible causes of incorrect landslide timings 280
In all of our case studies, our methods assign the wrong dates to small number of landslides. There are several possible reasons for this. There may be real changes in the time series that are not landslides, for example snowfall or melt, change in vegetation, change in soil moisture or human activity, which may be related to the landslide, for example the removal of material from a blocked road. Random noise in the SAR signal may also result in false landslide timings. We note that for future applications, the timing confidence within a landslide population can be separated into landslides timed by 3 or more methods and those 285 timed by only 2 methods (Table 2) Another possibility is that delayed or multi-stage failure occurred for some landslides. Our methods are designed to detect only a single failure. In the case where multi-stage failure results in more than one step change in the time series, the convolution in Sect. 2.5 will detect only the largest step change. Though it is beyond the scope of this study, it is in theory possible to assess if the time-series contain a second peak, of similar magnitude to the largest one, to assess the likelihood of multistage failure 290 or landslide reactivation.
Delayed failure seems particularly likely for Zimbabwe and Hiroshima, where a large proportion of the incorrect landslide timings are made up of the date pair immediately after the rainfall event (Fig. 4d, e, h). It is possible that some of the landslides in these inventories did not fail immediately during the rainfall, but instead failed after a delay of a few days due to rising pore pressure following rainfall infiltration within the hillslope (Iverson, 2000). This is particularly possible in the case of Z079D, 295 where the end of the rainfall event on 19 March 2019 coincides with the acquisition of the first post-event image, so that only a short delay would be required for the landslide to occur during the time window immediately after the rainfall (19-31 March 2019) rather than during the time window that spans the rainfall (7-19 March 2019). If these landslides are counted as correct in our analysis, the combined success rate in Zimbabwe is increased from 73% to 81%, bringing it in line with Hiroshima and Trishuli (Table 2), while for landslides timed by 3 or more methods ("⩾3M" in Table 2), the success rate is increased from 300 76% to 94%.
Although the Gorkha earthquake was followed by a large aftershock (12 May) and by the monsoon (approximate onset 9 June) (Williams et al., 2018), we are more confident of the true date of the landslides for this event. It is possible that some landslides could have been either triggered or reactivated by monsoon rainfall. However, none of the incorrect landslide timings in Nepal are in June, making this unlikely (Fig. 4a-c,f).

Why do some landslides have no timing estimation?
In all our case studies, a large proportion of landslides are not assigned any date by the amplitude methods. Some of these landslides, primarily in Nepal, lie in areas of foreshortening or layover in the SAR images and so were removed from the analysis (Sect. 2.3). Beyond this, there are several reasons that landslides may not be assigned a timing.
In Sect. 2.5, we discarded predictions for which the peak was too small in our convolution function. This improved the 310 specificity of our methods, but also required some correct predictions to be discarded. The reason for a small peak in the convolution function is either that the step change in the metric was too small, or the rest of the time series was too noisy. In Sect. 3.3, we demonstrated that landslides are less likely to be assigned a date if they occur in less vegetated areas or on slopes unfavourably oriented towards the sensor. This means that using SAR to detect landslides in arid areas remains a challenge, although methods based on coherence time-series may be more appropriate in such cases (e.g. Cabré et al., 2020). In the case 315 of the Geometric Shadows Method, there is no guarantee that a landslide will contain any shadow pixels, in which case no date can be detected using this method. This was observed to be more likely on slopes facing away from the sensor (Fig. 6).
Landslides will also have no signal in our methods if they are not located correctly in the SAR image, either due to distortion of the radar image or to inaccuracies in the optically-derived inventory. Milledge et al. (2021) and Pokharel et al. (2021) both observed differences in landslide shape and location when comparing different optically-derived inventories of landslides 320 triggered by the Gorkha earthquake, including the inventory of Roback et al. (2018) used here. Milledge et al. (2021) observed low spatial agreement between different inventories for all five of the trigger events they examined, suggesting this effect to be relatively common. Such differences can be due to problems georeferencing the optical imagery, imagery from different sources being used to compile different inventories and different teams of people carrying out the mapping. The landslide timing methods we present here require precise landslide locations. For the Geometric Shadows Method, which we expected 325 to be particularly sensitive to this form of spatial disagreement, we expanded the area covered by each landslide polygon by 10 m in an attempt to address this. However any mismatch beyond this scale between the landslide inventory and the SAR images is likely to lead to landslides not being assigned a timing.

The effect of shortening the time window
We used a "co-event" window of 6 months when testing the landslide timing methods in this paper. This time period was 330 selected to be roughly the duration of the Nepal monsoon. However, some applications, for example the case of successive storms, may not require such a long window. It is therefore useful to assess how the length of this time window affects the accuracy of predictions. In order to assess this, we took the tracks with the most complete time series (Z072A, H090D, Tr019D and BG019D) and assessed their performance over 2-8 month periods. Fig. 7 shows the percentage of assigned timings that are correct based on at least 2 methods for each track at each time period.

335
On all three tracks, particularly BG019D and H090D, the accuracy decreased as the co-event period was decreased. This was especially observed for periods of less than 5 months. We suggest that noise may be less attenuated in a shorter time series, resulting in increased numbers of false positives. Therefore, when applying our methods to future events, we recommend keeping a time window of at least 5 months even in the case where the landslide timing is known more precisely than this beforehand. This may explain the relatively poor performance in Bhote Kosi compared to the other case study areas, since 340 comparatively few images were available for this case study (Fig. 4). The loss of accuracy is recovered when the co-event period is further decreased to 2 months, possibly due to the comparatively small number of possible wrong date pairs available within a 2-month period.

InSAR Coherence
Interferometric SAR (InSAR) coherence is a measure of the signal quality of an interferogram (an image used to measure 345 ground deformation formed from two SAR images acquired over the same area at different times). InSAR coherence is sensitive to changes at the ground surface between the acquisition of the two SAR images: areas where the scatterers have changed significantly have high levels of noise in an interferogram and so a low coherence. Coherence is therefore sensitive to landslides and has previously been used to detect landslide densities or individual large landslides (Burrows et al., 2019(Burrows et al., , 2020Goorabi, 2020;Yun et al., 2015).

350
The coherence of each pixel in an interferogram can be estimated from the similarity in amplitude and phase change between the two SAR images for small groups of neighbouring pixels. Coherence surfaces and the phase data required for their calculation are not available through Google Earth Engine. However, coherence for Track 19 in Nepal has previously been calculated by Burrows et al. (2019), covering the Buri Gandaki and Trishuli inventories tested here. This allows us to compare methods of landslide timing based on SAR amplitude and InSAR coherence for these two case studies.
355 Burrows et al. (2019) processed the data at the same resolution used here (20 × 20 m) and used a 3 × 3 moving window to estimate coherence, so that the coherence surface has a resolution of 60 × 60 m. Similarly to the Landslide-Background Method, we obtained the median coherence of pixels within each landslide through time and the median coherence of pixels . within a 60-500 m buffer of each landslide polygon to give a background coherence. We then examined the ratio between the landslide and background coherence through time. Using this ratio performed better than using the landslide coherence 360 alone, probably because other factors, such as the length of time between the two images used to form the interferogram, can also effect coherence. Fig. 8 shows the median coherence ratio of a single landslide for different image pairs through time. This demonstrates two effects that we expect to see. First, the coherence that spans the landslide timing is low. This drop in coherence has previously been used to detect landslide locations (Burrows et al., 2019;Goorabi, 2020;Yun et al., 2015). However Sentinel-1 often has a low background coherence in vegetated areas due to its wavelength, which can make 365 any coherence decrease due to a landslide difficult to detect. Second, the coherence of post-event image pairs is higher than pre-event image pairs due to the removal of vegetation by the landslide (previously used by Burrows et al., 2020). Based on these two observations, we propose two landslide timing detection methods based on InSAR coherence time series.
Method C1: A step increase in the coherence ratio corresponds to the first post-event image pair.
Method C2: A temporary decrease in the coherence ratio corresponds to the co-event image pair. For each coherence pair, 370 this temporary decrease is calculated from the sum of the decrease in coherence ratio from the previous image pair to this one and the increase in coherence ratio from this image to the next (adapted from the ∆C_sum method of Burrows et al., 2020).
Overall, the coherence-based methods have a lower success rate than the amplitude-based methods (Table 2), indicating that incorporating these data would decrease the specificity of our methods. However, it is worth noting that of the 47 landslides correctly timed across the two events using the C1 and C2 combined, only 3 had already been timed using the combined 375 amplitude-based methods in Sect. 3.1, suggesting that the incorporation of coherence methods could increase sensitivity, if these could be made more reliable. Currently, only the Sentinel-1 SAR constellation acquires SAR data with sufficient coverage and acquisition frequency for use in landslide timing studies. These data are acquired at C-band, which usually has low coherence in vegetated areas. Lband data is better suited to InSAR-coherence-based landslide detection in vegetated areas (Burrows et al., 2020). The planned 380 NASA-ISRO NISAR mission has a similar acquisition strategy to Sentinel and will acquire L-band SAR. It will be worth reassessing the potential of InSAR coherence time series for landslide timing detection following the launch of this satellite.

Application to future events
We have developed methods that allow around 20% of landslides in an inventory to be dated at the 6-12 days scale with ∼ 80% confidence. In the case of multiple successive storms, this timescale should be sufficient to attribute landslides to a given 385 event. For monsoon landslide dating, this timing is not sufficiently precise to allow building of intensity-duration or intensityantecedent rainfall thresholds (e.g. Bogaard and Greco, 2018). However, it should allow us to establish whether landslides occur in temporal clusters that relate to specific peaks in rainfall or are distributed throughout the monsoon. These two endmembers would have very different implications in terms of hydrological and slope stability modeling, and thus on hazard evaluation. It should also allow us to better constrain whether landslides systematically occur with a specific delay after the 390 onset of the monsoon and/or simultaneously with reported flooding or bursts of intense rainfall (Gabet et al., 2004). In both cases, since relatively few landslides are timed by SAR methods, some inference will be needed to attribute time stamps to neighbouring landslides lacking a SAR signal. The SAR methods could also be combined with other methods of establishing landslide timing if these are available e.g. small gaps in cloud in optical satellite imagery, reports of individual landslides or seismic data (Bell et al., 2021;Hibert et al., 2019;Yamada et al., 2012).

395
An algorithm could be designed to automatically detect both the location and timings of rainfall-triggered landslides by combining SAR-based timing methods with automated landslide mapping methods based on multi-spectral imagery (e.g. Behling et al., 2014;Ghorbanzadeh et al., 2021;Jelének and Kopačková-Strnadová, 2021;Milledge et al., 2021). They could also be applied to other forms of landslide inventory, for example those based on LiDAR scans or high resolution optical images that allow landslide volumes to be estimated.
In the case of prolonged rainfall events, landslide inventories compiled from optical satellite imagery are often poorly constrained in time. Here we present methods of using Sentinel-1 SAR images in Google Earth Engine to identify the timing of rainfall triggered landslides to within a few days. We find that by combining ascending and descending track SAR, it is possible to date ∼ 20% of landslides in an inventory with an accuracy of ∼ 80%. A small number of landslides (3-8%) can be timed 405 with accuracy of 90-95%. These methods will allow us to generate multi-temporal landslide inventories for long rainfall events, unlocking comparisons between rainfall data, hydrological models and triggered landsliding. It also suggests that SAR amplitude time-series could be combined with multispectral imagery in an algorithm aimed at detecting and mapping landsliding over regional scales.
Code and data availability. Sentinel-1 GRD and Sentinel-2 data are available open-access from ESA Copernicus and were accessed through evolution of rainfall-triggered landslides using radar and optical satellite data". We thank Robert Emberson for sharing some landslide inventories.