A data-driven evaluation of post-fire landslide susceptibility

Wildfires change the hydrologic and geomorphic response of watersheds, which has been associated with cascading hazards that include shallow landslides and debris flows. This study evaluates post-wildfire landslide trigger characteristics by comparing precipitation preceding landslides at both burned and unburned locations. Landslide events are selected from the NASA Global Landslide Catalog (GLC) to facilitate regional inter-comparison. Fire and precipitation histories for each site are established using MODIS global burned area and CHIRPS precipitation data, respectively. Analysis of normalized 5 seven-day accumulated precipitation for sites across all regions shows that, globally, landslides at burned sites are preceded by less precipitation than landslides without antecedent burn events. This supports the hypothesis that fire increases rainfalldriven landslide hazards. An analysis of the seasonality of landslides at burned and unburned locations shows that landslidetriggering storms in burned locations tend to exhibit different seasonality from other rainfall-triggered landslides, with a variety of seasonal shifts ranging from approximately six months in the Pacific Northwest of North America to one week in the 10 Himalaya region. Overall, this manuscript offers an exploration of regional differences in the characteristics of rainfall-triggered landslides over a broad spatial scale and encompassing a variety of climates, geographies, and burn conditions.


Sources and methods for landslide data collection
It is resource-prohibitive to conduct a continuous systematic search for mass movements either in the field or with satellite 90 observations. As a result, many of the most accurate and complete methods for systematically identifying landslides can presently only be used over limited spatial and temporal domains. For example, Lee and Pradhan (2007) identified landslides from aerial photograph interpretation and a field survey over the ∼ 800 km 2 Selangor area in Malaysia, and Nefeslioglu et al. (2010), used an inventory based on aerial photographs taken in 1955-1956 to analyze landslide susceptibility over a ∼ 175 km 2 area near Istanbul, Turkey. An alternative to manual identification either in the field or using photographs is automatic or semi-95 automatic landslide detection using image processing on aerial imagery, LiDAR surveys, or Synthetic Aperture Radar (SAR).
These automated methods are typically applied over similarly small domains due to challenges with obtaining imagery and compiling training datasets. For example, Martha et al. (2013) used aerial imagery over ∼ 120 km 2 in the Himalayas, while Mezaal et al. (2017) used LiDAR over the 26.7 km 2 Cameron Highlands of Malaysia. SAR interferometry can be used for identification of pre-landslide motion, as was done by Lu et al. (2012) over the ∼ 1500 km 2 Arno basin in Italy. In addition, 100 several SAR techniques have been employed to identify post-landslide scars, including SAR amplitude mapping of landslides triggered by the Gorkha, Nepal earthquake in 2015 over a 14, 500 km 2 area (Meena and Tavakkoli Piralilou, 2019), coherence mapping of interferometric SAR the same earthquake-triggered landslides (Burrows et al., 2019), and the wildfire-triggered landslides over ∼ 60 km 2 of the area burned by the 2017 Thomas Fire in California (Donnellan et al., 2018). While automated landslide detection as deployed in the above studies is continually undergoing promising advances, at the time of this analysis 105 it has not yet been used to compile an inventory over a broad enough spatial domain to facilitate inter-regional comparisons.
Such records collected in an uncoordinated effort over small domains are unsuitable for regional inter-comparisons such as we have undertaken here because these records do not contain standardized information for every region, are often unpublished (van Westen et al., 2006), and are unlikely to have daily temporal resolution that would allow comparison with the precipitation record (Kirschbaum et al., 2010). 110 For this study, we chose to use the NASA Global Landslide Catelog (GLC, Kirschbaum et al., 2010). As with the few other regional and global databases available, the broad domain of the GLC comes coupled with issues of precision and completeness. The sources of GLC data are second-hand observations made by the news media, governmental organizations such as departments of transportation, and some available scientific reports (Kirschbaum et al., 2010). The absence of a systematic search for landslides across the entire database domain results in a substantial spatial bias towards populated areas where 115 landslides happen to be noticeable. News reports also suffer relatively high location uncertainty (as much as 50 km) depending on how specific the source article is about the location (Kirschbaum et al., 2010). Finally, though the GLC does contain some information about the landslide mechanisms that would allow landslides to be classified, for example, as debris flows or shallow landslides, these data may not be reliable since they are not based on direct observations. Despite limitations in accuracy and completeness, the GLC was chosen for this study primarily because as of this writing it offers the largest spatial 120 and temporal range of any catalog. The GLC contains a large sample of landslides at unburned control basins from across the globe (n = 5313) as well as a substantial proportion of landslides identified in this study as having occurred in recently burned areas (n = 489; 9.2 %).

Towards a global picture of global landslides susceptibility
This study seeks to test the hypothesis that wildfires increase landslide susceptibility by evaluating antecedent precipitation 125 at both burned and unburned landslides locations. Some existing local and regional studies (Cannon et al., 2010;Rupert et al., 2003) have assessed the impact of wildfire on landslide susceptibility, but have not included unburned locations in their analyses. Other studies have also featured the GLC data and a global spatial extent, with a focus on validating large-scale landslide hazard models (Kirschbaum and Stanley, 2018). This analysis is unique from other regional and global studies in that it combines the broad scope of the GLC data with an exploration of the role of wildfire in landslide susceptibility. This 130 study is also distinct from others that focus on the role of wildfire on landslide sites (Gartner et al., 2009) in that here burned sites are contrasted with unburned sites instead of previous observations of the same location. Finally, in contrast to postwildfire landslide studies focused on a specific regions like the western US (Cannon and DeGraff, 2009), southern California (Gartner et al., 2014), or southeast Australia (Nyman et al., 2011), this study combines the GLC with globally-observed fire and precipitation data to offer unique insights into the role of fire on landslide susceptibility in diverse regions across the globe.

Methods
We first describe the landslide data (Sect. 2.1), the study regions (Sect. 2.2) and fire data (Sect. 2.3). The precipitation data (Sect. 2.4) leading up to the date of each landslide were compared using three approaches. First, the seven-day running total precipitation depth percentile (see Sect. 2.4) in the weeks preceding each landslide was used as a proxy for landslide susceptibility. This percentile value was compared between burned and unburned sites within each region and for all included landslides 140 (see Sect. 2.5). Next, seven-day precipitation percentiles were compared with bootstrapped samples from burned and unburned sites separately (see Sect. 2.6) to confirm the findings from the depth percentile analysis and also to draw out differences in storm timing between burned and unburned groups. Finally, the precipitation frequency in the burned and unburned groups in the months and years surrounding each landslide (see Sect. 2.7) was examined to identify shifts in the seasonality of landslides at burned sites relative to the unburned group. These seasonality results were augmented with kernel density estimates 145 of landslide occurrence by day-of-year at burned and unburned sites for each region.

Landslide data
A large sample (n = 5313) of rainfall-triggered landslides was obtained from the GLC. Landslide locations are shown in Fig.   1, along with a summary of fire and precipitation information obtained for those locations from the sources listed in Table   1 (see Sects. 2.3 and 2.4). The GLC provides a large collection of events taking place in a variety of climates such that, in 150 combination with spatially continuous observations of fire (500m Moderate Resolution Imaging Spectroradiometer [MODIS] Burned Area by Giglio et al., 2018) and precipitation (5.5km Climate Hazards group InfraRed Precipitation with Station data burned and unburned groups in the regional insets; burned/unburned classification at the time of the landslide in (c) and (d), with regional insets showing kernel density portrayal of the fraction of burned area for the landslide locations from the three years preceding the landslide; and the precipitation percentile on the day of the landslide in (e) and (f), with regional insets of kernel density estimates (violin plots) of the climatological (1981 − 2020) seasonal precipitation magnitude (mm) including a reference line indicating the median seasonal average across all sites globally. Country boundaries were obtained from the maps R package (Deckmyn et al., 2018) 6 https://doi.org/10.5194/nhess-2021-111 Preprint. Discussion started: 28 April 2021 c Author(s) 2021. CC BY 4.0 License.
[CHIRPS] by Funk et al., 2015) data, it is well suited for comparing the diverse precursors under which post-wildfire landslides occur.
In order to reduce errors resulting from including a variety of types of rainfall-triggered landslides within the same dataset, 155 the selected landslides were limited to those categorized by a 'landslide trigger' value of 'rain,' 'downpour,' 'flooding,' or 'continuous rain.' Landslides with a second trigger such as an earthquake were eliminated. Snowmelt-driven landslides were also not included because the impact of precipitation is delayed in those cases -an analysis of the snow record in California/Nevada revealed only a single event with enough antecedent snow to suggest it could have been mislabeled. Only records that were indicated to have occurred within a radius of 10 km or less of the recorded location were included, since the landslides with 160 lower location accuracy presented problems for wildfire classification. In addition, only landslides between 50 • S and 50 • N latitude were included, and the events occurring before the year 2000 were omitted, so as to ensure coverage by both fire and precipitation datasets (see Table 1).

Study regions
To compare the differences in landslide triggers in different climates, we divided the landslides into regions (see Fig. 1 panels 165 (a)and (b)). Regions were determined using the AGglomerative NESting (AGNES) hierarchical clustering algorithm (Kaufman and Rousseeuw, 2009) considering the latitude and longitude of the landslides, and clusters were subsequently combined, split, or eliminated on the basis of sample sizes as described below. First, the cluster tree was truncated at 30 clusters, after which all the clusters with fewer than 100 data points or less than 5 % burned sites were eliminated. Cases where two nearby regions with lower numbers of landslides, for example, Central America and Caribbean/Venezuela, were joined manually. Finally, 170 the largest region, encompassing Western US and Canada, was split into three sub-regions based on an additional identical clustering process over this sub-domain. The final regions are shown in Fig. 1

panel (a). The Pacific Northwest of North
America was included even though the percentage of burned sites is lower than threshold, but at 4.4 % it was nearly double the highest percentage among the eliminated regions (2.25 % in the Eastern US). Some landslides were not included in any of the final regions. These events were not, however, eliminated from any analysis of all landslides.

Fire data
For each landslide, a circle centered at the landslide location and with a radius of the location accuracy was computed and each 500 m pixel within the circle was extracted from the MODIS Burned Area dataset (Giglio et al., 2018). Fire affects the landscape over a large range of temporal scales in different settings. Previous studies suggest that the post-wildfire increase in landslide susceptibility peaks within the first six months, but that a second time period of increased susceptibility can appear 180 at 3 years or even longer (DeGraff et al., 2015;Gartner et al., 2014). Landslides were classified as burned if any part of the area where the landslide occurred was burned at some point within the three years prior to the event to capture both waves of increased susceptibility without over-identifying landslides areas where fires occur every few years. The fraction of pixels that were burned over the 3-year antecedent period was then computed, and landslides classified as burned if there was any overlap

Precipitation data
Time series of precipitation at the landslide sites were obtained from the CHIRPS precipitation dataset (Funk et al., 2015).
Precipitation was averaged for each landslide location within the radius of the provided location accuracy. Additional preprocessing steps described below were performed to distinguish anomalously high precipitation events from potential seasonal shifts and climatic differences across sites.

200
A 7-day running average of antecedent precipitation was computed to enable direct comparison of landslides triggered by storms of different lengths and intensities. While including an estimate of the soil moisture was outside the scope of this study, 7-day antecedent rainfall indices have been used by other modelling studies as a surrogate for soil moisture in a combined indicator of landslide susceptibility (James and Roulet, 2009;Kirschbaum and Stanley, 2018). In addition to allowing a more equal comparison of landslide triggers which fall within throughout this spectrum of storm intensity, the running average 205 is less sensitive to small errors in precipitation and landslide date accuracy. Figure 1 panels (e) and (f) show these 7-day cumulative-precipitation percentiles, as well as the climatological seasonal average precipitation, revealing that the Western US is dominated by dry summers, while the lower-latitude regions exhibit wetter summers and in some cases monsoons.
We further screened landslide events with no recorded precipitation in advance of reportedly rainfall-triggered landslides. Figure 2 shows a quality control sub-analysis for the California/Nevada area to investigate the need for data screening on the 210 basis of inconsistencies between the landslide and precipitation record. This region was chosen for the quality control analysis, because of its high density of precipitation data and variety of climate conditions, useful for identifying erroneous landslide precipitation. We found 14 % (73 of 533) of the landslides in this region had no triggering precipitation event recorded in the CHIRPS data. Since the GLC contains only rainfall-triggered landslides, the lack of precipitation in these cases was likely a result of errors in either the precipitation data or landslide data. A comparison with the Daymet precipitation dataset over the 215 same domain revealed that the two precipitation datasets frequently did not agree on these zero-precipitation landslide events, suggesting that the problem largely originated from the precipitation data themselves. Furthermore, the concentration of data points on the x and y axes of Fig. 2 suggests that disagreements on precipitation occurrence are distinct from disagreements on the non-zero amounts of precipitation and potentially a separate source of error. To limit the effect of these inconsistent data points on the results, all landslides worldwide with no measured precipitation in the six days before and one day after the event 220 were removed from the global study (367 of 5680 or 6.5 % removed for a final n = 5313).

Precipitation percentile experiment 230
This experiment compares the 7-day precipitation percentile in the burned and unburned groups in the time leading up to a landslide. The percentile indicates the degree to which landslide-triggering storms were exceptionally large and also serves as a proxy for relative landslide susceptibility. A one-sided Mann-Whitney hypothesis test was used to ascertain whether the precipitation percentiles of burned sites were less than the precipitation percentiles of unburned sites. Deviations between the burned and unburned groups defined by a p-value less than 0.05 on the Mann-Whitney test indicate statistically significant 235 differences in the landslide susceptibility of the two groups. The null hypothesis of the Mann-Whitney test was that the median percentile of the burned sites is greater than or equal to the median percentile of the unburned sites. The precipitation per-centile distributions of all 7-day periods with non-zero precipitation are uniformly distributed as a result of the pre-processing described in Sect. 2.4, making the Mann-Whitney test the most appropriate hypothesis test. However, since zero-precipitation periods are excluded, this method cannot account for differences in the frequency of precipitation across different climates, but 240 rather reflects differences in the magnitude of 7-day precipitation totals.

Bootstrapped samples experiment
In order to evaluate how anomalous the precipitation events preceding burned and unburned landslides were to "typical" local climate conditions at the landslide locations, we compared them to bootstrapped samples from other years to obtain a clearer signal. One hundred samples were taken from the 38-year precipitation records to match the locations and DOY of the observed 245 landslides, but from randomly selected years (n ≈ 100 for the individual samples, with the actual sample number adjusted by region so that all sites were selected evenly). Sampling was repeated for burned and unburned groups within each region as well as for all the landslides in the study. These samples are representative of precipitation for a particular lead time at the landslide locations and serve as a control dataset with which to compare the pre-landslide precipitation. Next, the observed landslide-year precipitation across all sites in the group was tested against each bootstrap sample using a Mann-Whitney test, 250 with the null hypothesis that the sample median was less than or equal to the median of the precipitation percentiles from that day of the year in the entire record from 1981-2020. This produced a distribution of p-values that represent the likelihood that the precipitation leading up to the landslides varied from the control baseline.
This sampling method, though more complex, helps to reduce noise in the hypothesis test results due to different sample sizes in different regions. It also provides more information on general landslide susceptibility of each region rather than only 255 the relative susceptibility of burned and unburned sites. Finally, it includes measurements of zero precipitation, which were eliminated from the direct comparison because of long-term climatic differences in precipitation frequency between burned and unburned sites in all regions.

Landslide seasonality experiment
To investigate the potential role of wildfires in affecting landslide seasonality, we estimated precipitation frequency at the 260 landslide sites over time by computing the fraction of sites in the burned and unburned groups that had precipitation on any given day. As with the percentiles and the bootstrap p-values, frequency estimates were computed relative to the landslide event rather than by calendar date, resulting in time coordinates measured in 'years before the landslide'. Precipitation frequency was estimated for two years before and after the landslide in order to highlight changes in the magnitude and phase of the precipitation pattern. We found that in most regions there was a long-term difference in the mean annual frequency, likely 265 because fires occur more often in areas with drier climates (Liu et al., 2014). These persistent differences between burned and unburned sites were removed by subtracting the mean precipitation frequency for both the burned and unburned groups. Finally, we took a 90-day running average to reduce noise in the data and thereby make it easier to visually identify any long-term shifts in landslide occurrence. These frequency estimates are not normalized by season, which means that unlike the previous two to annual precipitation cycles.
Additional seasonality analysis was performed to provide insight into the times of year that landslides occur at burned versus unburned sites. Kernel density estimates of landslide occurrence throughout the year were compared between the burned and unburned groups. This seasonality analysis would highlight a shift from Fall to Spring but, in contrast with the frequency analysis, it does not indicate the precipitation conditions under which landslides typically occur. Together, the frequency and 275 seasonality analyses can show both the seasonal shift as well as any changes in landslide occurrence relative to annual precipitation patterns.

Precipitation percentile experiment
The distributions of precipitation event percentiles for all the included landslides are shown in Fig. 3. The precipitation per-280 centile increases for all groups as the date of the landslide approaches, confirming that these rainfall-triggered landslides are generally preceded by an increase in total precipitation depth. Notably, when considering all landslides together (Fig. 3) the precipitation events that triggered landslides at burned sites were significantly smaller than those that triggered landslides at unburned locations (Mann-Whitney test, 95 % confidence). This difference supports the overarching hypothesis that wildfire does in fact increase landslide susceptibility, since landslides in the period after a fire can be triggered by less precipitation than 285 might normally be required to cause mass movement. An examination of each region separately reveals that the difference in precipitation percentile between burned and unburned sites is present in some regions but not in others (see Fig. 3). The California area ( Fig. 3 panel (b)) has a particularly strong signal, whereas tropical regions do not show any significant differences between precipitation at burned and unburned sites or display the reverse effect of higher precipitation percentiles for unburned locations than burned locations. orange or purple curve was tested against the black curve to obtain the boxplots of p-values at 0 days before the landslide.
A clear difference between burned and unburned sites is shown for the same regions as in Fig. 3, but with the addition of Southeast Asia. Beyond the emergence of a signal in Southeast Asia, additional differences between regions in the timing of  samples were drawn from the same DOY and locations as the landslides but from a randomly selected year. In panels (a)-(g), box plots of p-values represent the degree to which the landslide-triggering precipitation differed from climatological precipitation with lower values indicating a larger difference between the two precipitation distributions. To illustrate whether there are differences between burned and unburned groups when the p-values below 0.05 (significant at 95 % confidence, shown as a dashed black line), the y-axes are shown with a probit transform to expand that section of the axis. The y-axis has also been inverted so that larger differences in precipitation (lower p-values) are higher on the y-axis for consistency with the percentile plots in Fig. 3. In panels (h)-(u), an example of the kernel density estimate (kde) for day-of-landslide precipitation in black separated by burned and unburned groups is compared with kdes of all bootstrapped samples in orange (burned group) or purple (unburned group). percentile analysis.
Different storm timing is apparent among the regions, with implications for potentially region-specific physical processes associated with landslide triggers. Firstly, in the Himalayas and Southeast Asia (Fig. 4 panels (f) and (g)) precipitation rises at a similar rate for each group, indicating that landslides at burned and unburned locations are triggered by similar precipitation increases. Curiously, the bootstrap analysis reveals a long-term difference between burned sites and unburned sites in the 305 Mann-Whitney p-value for Southeast Asia despite location-specific normalization, suggesting that the landslides at unburned locations might be primarily triggered in years that are wetter than usual on a monthly or possibly seasonal scale. In the Pacific Northwest and California ( Fig. 4 panels (d) and (b)), the burned sites exhibit shorter but more intense storms than the unburned sites in the week preceding the landslide. Under the assumption that shorter, more intense storms are associated with runoffdriven landslides while longer storms that allow more time for the soil column to saturate are associated with infiltration-driven 310 landslides, this difference in storm timing suggests that the burned landslide locations are largely runoff-driven while rainfalltriggered landslides at unburned locations are infiltration-driven (Cannon and Gartner, 2005). The Mann-Whitney p-values for the burned group remain well above 0.05 just days before the landslide as the p-value for the unburned group begins to fall.
In the Intermountain West ( Fig. 4 panel (c)) antecedent precipitation for the burned group is generally characterized by a dry spell going back thirty days or more. Twenty to thirty days before the landslide p-values are consistently above 0.9, suggesting 315 a high likelihood (> 90 %) that there was less precipitation than usual during that time. occur nearly year-round when considering all regions together, but the other panels in Fig. 5 show that within any particular region, fires occur only during a distinct time of year. However, the delay between fire and landslide is not consistently equal to the length of time between fire season and the following landslide season. The landslides are distributed such that 4854 % occur within a one year after the fire and at least 10 % of sites in each region the fire occurred between two and three years before the 325 landslide. Since both landslides and fires have seasonal patterns, the typical delay between fire and landslide for each region appears to be primarily related to the relationship between fire season and landslide season. For example, California has a long fire season and a shorter landslide season, and so when fires occur at the end of winter, immediately after landslide season, there is typically a longer delay before the landslide than when fires occur immediately before landslide season. By contrast, in the Himalayas the delay between fire and landslide is relatively uniform due to a shorter fire season that does not overlap with 330 the landslide season. Figure 6 shows differences in seasonality between burned and unburned landslide seasonality on the right and the results of the precipitation frequency analysis on the left. While all regions except for Central America (Fig. 6 panel (l)) display some