Spatial distributions of earthquake-induced landslides and hillslope preconditioning in the northwest South Island, New Zealand

. Current models to explain regional-scale landslide events are not able to account for the possible effects of the legacy of previous earthquakes, which have triggered landslides in the past and are known to drive damage accumulation in brittle hillslope materials. This paper tests the hypothesis that spatial distributions of earthquake-induced landslides are determined by both the conditions at the time of the triggering earthquake (time-independent factors) and the legacy of past events (time-dependent factors). To explore this, we undertake an analysis of failures triggered by the 1929 Buller and 1968 Inangahua earthquakes, in the northwest South Island of New Zealand. The spatial extents of landslides triggered by these events were in part coincident. Spatial distributions of earthquake-triggered landslides are determined by a combination of earthquake and local characteristics, which inﬂuence the dynamic response of hillslopes. To identify the inﬂuence of a legacy from past events, we ﬁrst use logistic regression to control for the effects of time-independent variables. Through this analysis we ﬁnd that seismic ground motion, hillslope gradient, lithology, and the effects of topographic ampliﬁcation caused by ridge-and slope-scale topography exhibit a consistent inﬂuence on the spatial distribution of landslides in both earthquakes. We then assess whether variability unexplained by these variables may be attributed to the legacy of past events. Our results suggest that hillslopes in regions that experienced strong ground motions in 1929 were more likely to fail in 1968 than would be expected on the basis of time-independent factors alone. This effect is consistent with our hypothesis that unfailed hillslopes in the 1929 earthquake were weakened by damage accumulated during this earthquake and its associated aftershock sequence, which inﬂu-enced the behaviour of the landscape in the 1968 earthquake. While our results are tentative, they suggest that the damage legacy of large earthquakes may persist in parts of the landscape for much longer than observed sub-decadal periods of post-seismic landslide activity and sediment evacuation. Consequently, a lack of knowledge of the damage state of hillslopes in a landscape potentially represents an important source of uncertainty when assessing landslide susceptibility. Constraining the damage history of hillslopes, through analysis of historical events, therefore provides a potential means of reducing this uncertainty.


Introduction
Regional landslide-hazard assessments rely on models that upscale our conceptual understanding of fundamental controls on landslides, through analysis of the influence of different proxy variables on landslide occurrence (e.g.Capolongo et al., 2002;Garcia-Rodriguez et al., 2008).Most studies to date have addressed spatial correlations between the distribution of landslides and variables that provide proxies for seismic ground motions and the modelled stability of hillslopes (e.g.Dai et al., 2011;Meunier et al., 2007Meunier et al., , 2008;;Kritikos et al., 2015).These studies implicitly rely upon a static model of hillslope sensitivity to landslide triggering.In other words, the predicted number of landslides triggered by any given trigger event, or the susceptibility to landsliding in that event, will not vary through time.However, this assumption is at odds with observations of increased rainfall-triggered landslide activity above baseline rates observed in the wake of large earthquakes (Hovius et al., 2011;Saba et al., 2010;Tang et al., 2011;Dadson et al., 2004).Similarly, data from the 2010-2011 Canterbury earthquake sequence reveal landslide triggering at lower ground accelerations following the February 2011 earthquake, which caused cracks to develop in hillslopes that subsequently failed in later earthquakes in the sequence (Massey et al., 2014a, b;Mcfadden et al., 2005).These observations suggest that hillslopes may retain damage from past earthquakes, which makes them more susceptible to failure in future triggering events.Note that here we define failure as the total collapse of a hillslope where the failed mass evacuates the failure plane and moves downslope to leave a discernable, bare-earth scar.According to the classification of Keefer (1984Keefer ( , 2002)), these types of failures are generally grouped as disrupted slides, given the significant internal disruption exhibited by the landslide mass.In the classification system of Varnes (1978;updated by Hungr et al., 2014), this group includes rock and debris falls and slides, as well as rock avalanches.Globally, disrupted landslides are estimated to comprise the majority, ∼ 86 %, of reported earthquake-induced landslides (Keefer, 1984(Keefer, , 2002;;Mcfadden et al., 2005).
One mechanism by which hillslopes could be made more susceptible to failure is progressive brittle damage accumulation in hillslope materials, whereby permanent slope displacement leads to cracking and dilation of the mass (Petley et al., 2005;Nara et al., 2011;Bagde and Petroš, 2009;Li et al., 1992b).Damage accumulation occurs near the surface within hillslopes (Clarke and Burbank, 2011), as gravitational stress coupled with seismically and hydrologically induced changes in the stress distribution cyclically load and unload hillslope materials.Cyclic loading at stresses lower than the static strength of the material can produce irreversible localised strain damage (Suresh, 1998;Schijve, 2001), in natural rocks with abundant pre-existing microcracks or pores (Attewell and Farmer, 1973;Li et al., 1992a;Bagde and Petroš, 2005).While the effect of an individual cycle may be insignificant, the nonlinear accumulation of damage over repeated cycles can lead to the eventual failure of the material, via a progressive mechanism of failure (Petley et al., 2005;Leroueil et al., 2012).Brittle deformation of this type has been observed in soil at low confining pressures (1-250 kPa), in mudrocks at confining pressures up to 2 MPa, and at greater confining pressures in harder geological materials (Petley and Allison, 1997;Evans et al., 2013).As this mechanism occurs in the fabric of brittle rock or cohesive soils (bonded or cemented materials, where strain is localised during failure), it is likely to be common to most disrupted types of landslide induced by earthquakes.Exceptions to this are shallow colluvial failures in cohesionless soil and cases of failure in very poor quality, soft rock masses or soft layers, where material behaves in a ductile manner (Petley and Allison, 1997).Where earthquake-induced landslide failure develops progressively, via brittle deformation, hillslopes may retain damage from past earthquakes.Whether or not a hillslope fails in response to an earthquake will be a function of both the current event and, by definition, the history of damage accumulated in that hillslope from previous events.The absence of this historical information from landslide analyses and predictive models potentially represents a significant gap in our understanding of factors that control the distribution of landsliding.
If damage from previous earthquakes does influence patterns of landsliding in subsequent earthquakes, then it is reasonable to hypothesise that spatial distributions of landslides should be at least partially correlated with the ground motions from past earthquakes.In order to investigate the role of hillslope damage history in conditioning landslide distributions, we test this hypothesis through analysis of the spatial distribution of landslides triggered by two large (M w > 7) earthquakes which occurred in close proximity in the northwest South Island of New Zealand.First, we present inventories of landslides triggered by the 1929 Buller and 1968 Inangahua earthquakes.Second, we undertake a spatial analysis of the distributions of both events, using logistic regression.Third, we use the results of this analysis to test the influence of the 1929 earthquake on the distribution of landslides triggered by the 1968 earthquake.

The 1929 Buller and 1968 Inangahua earthquakes
The 17 June 1929 Buller (Murchison) earthquake (M w = 7.7; Dowrick and Rhoades, 1998;Dowrick, 1994) and the 24 May 1968 Inangahua earthquake (M w = 7.1; Anderson et al., 1994) both triggered landslides over a large area (Fig. 1).The epicentres of the two earthquakes were ∼ 21 km apart, whilst at their closest point the mapped surface expressions of the coseismic faults lie 7 km apart.Source analyses described below indicate that the earthquakes have very similar reverse thrust focal mechanisms, with small components of left-lateral strike-slip (Dowrick and Smith, 1990 son et al., 1993, 1994;Stirling et al., 2007), isoseismal contours (Dowrick, 1994;Adams et al., 1968), mapping coverage regions and earthquake-induced landslides mapped in this study.et al., 1994).Isoseismal maps (Dowrick, 1994;Adams et al., 1968) suggest that ground motions from the two events had a MMI VIII overlap area of ∼ 3505 km 2 and a MMI IX overlap area of ∼ 584 km 2 (Fig. 1).

Coseismic sources and ground motion
The White Creek Fault has been identified as the source of the 1929 earthquake, although surface faulting was only observed along an 8 km length of the fault (Fyfe, 1929;Hender-son, 1937).Back analysis of seismic data (Doser et al., 1999), ground motion intensities (Dowrick, 1994), and coseismic landslides (Pearce and Oloughlin, 1985;Hancox et al., 2002) suggests a unilateral rupture extending 30-50 km to the north of the epicentre.This corresponds with the mapped geological (ground surface) trace of the White Creek Fault.Estimates of dip angle range from 60-70 • based on surface displacement observations (Henderson, 1937) to 46 • ± 13 • based on inversion of data from seismic stations (Doser et al., 1999), and 45 • based on elastic dislocation modelling (Haines, 1991).Doser et al. (1999) inferred a focal depth of 9 ± 3 km.To approximate the 1929 seismic source geometry in our analysis, we use the surface fault line and fault parameters of the White Creek Fault as used in the New Zealand probabilistic seismic hazard model (Stirling et al., 2000(Stirling et al., , 2002(Stirling et al., , 2007(Stirling et al., , 2012;;Berryman, 1980;Haines, 1991).This model assumes a fault plane striking 10 • , and dipping at 45 • from the surface to a maximum depth of 12 km, with a dip direction of 100 • .
The seismic source geometry of the 1968 earthquake has been constrained through an integrated geological, geodetic and seismological source model (Anderson et al., 1993(Anderson et al., , 1994)).We use a single fault plane trending northeast (25 • ), dipping at ∼ 45 • from a depth of 10-15 km to within ∼ 1 km of the surface (i.e.no primary ground surface rupture), with a dip direction of 295 • extending around 30 km in length (Anderson et al., 1993(Anderson et al., , 1994)).Earthquake parameters for both events are summarised in Table 1.
As coseismic landslide occurrence is driven by seismic shaking, it is important that we constrain the spatial pattern of ground accelerations.The strength of seismic ground accelerations attenuates with distance from the seismic source.However, the regional distribution of ground acceleration is also subject to the effect of rupture directivity and regional variation in the damping effect of earth materials (Abrahamson et al., 2008;Campbell and Bozorgnia, 2008).In an attempt to account for these effects in the case of the 1968 earthquake, we also make use of the USGS Shakemap output for this event (USGS, 2014).This Shakemap is based on the fault model described above and uses ground motion data from 15 seismic stations across New Zealand, 3 of which are within or just beyond the area of landslide mapping conducted here (Fig. 1), as well as estimates of peak ground acceleration (PGA) derived from reports at 159 additional sites.Although this model is still subject to uncertainty, by incor-  porating observed ground motions and site amplification factors, it can potentially provide a more accurate representation of the regional distribution of ground motion.PGA estimates derived from scratch-plate records at Reefton, Westport and Murchison report ground accelerations of 0.58, 0.30 and 0.36 g respectively (Adams et al., 1968;Dowrick and Sritharan, 1993), with which the Shakemap data set is consistent.

Earthquake-induced landslides
Both earthquakes triggered widespread landsliding throughout the area that experienced intensities of MMI = VIII to X.We review the types of landslides triggered by the earthquakes and outline our methodology for producing landslide inventories for the two events.

Landslide types
Most failures triggered by these earthquakes were disrupted rock and debris slides, rockfalls and rock avalanches, with very few coherent landslides and lateral spreads seen in the field or in aerial photos (Hancox et al., 2002(Hancox et al., , 2014)).In Figs.2-6 we present examples of these different landslide types from the two earthquakes.Note that an extended re- view of major landslides and landslide types is presented in Hancox et al. (2014).
Rockfalls were commonly triggered on steep scarps of Tertiary limestone, granite and greywacke, with numerous failures ranging from individual, small boulders to large falls of 10 5 m 3 (Fig. 2).Debris slides were the most frequent type of landslide triggered by the earthquakes and were common in areas of granite and greywacke (Fig. 3).Several examples of large rock avalanches were triggered by the earthquakes.The 1929 earthquake triggered the 18 million m 3 Lake Stanley rock avalanche (Fig. 4a), in Palaeozoic conglomerate and volcanics around 90 km north of the epicentre.Although this landslide is 35 km north of the present study area and is not included in the 1929 landslide data set (which covers only the southern half of the landslide-affected area) it is typical of the 10 largest landslides that occurred in 1929 (Hancox et al., 2002).The landslide is around ∼ 2 km long with an elevation range of 800 m.The largest landslide triggered by the 1968 earthquake was a 5 million m 3 rock avalanche (Fig. 4b).This failure occurred in weathered granite, running out about 1.2 km to the valley floor and about 100 m up the opposite side of the valley.
Several large rockslides were also triggered by the earthquakes.For example, the 1929 earthquake triggered the 18 million m 3 Matakitaki landslide (Hancox et al., 2002).This dip slope rockslide travelled ∼ 1 km across the valley floor, destroying two farm houses and killing four people, and formed a landslide dam (Fig. 5a). Figure 5b shows the intensity of landslide damage in the Matiri Valley, an area close to the seismogenic fault, where landslide scars from 1929 are still clearly visible today.The 1968 earthquake triggered the 3 million m 3 Oweka rockslide, a disrupted mass of muddy sandstone that fell from a vegetation-covered slope (Fig. 6a).The largest (2.8 million m 3 ) rotational landslide triggered by the 1968 earthquake occurred on a 100 m high terrace in sandy ("Blue Bottom") mudstone (Fig. 6b).
In most of these failure types, we might reasonably expect the process of material failure to involve some component of brittle deformation, given the low temperature and confining pressure in near-surface materials.Notable exceptions to this may include structurally controlled failures along ductile bedding planes.For example, field observations from the Oweka landslide suggest that, for a large semi-intact section of the landslide, the mechanism of movement was sliding on an extensive bedding plane coated with a thin layer of plastic clay (Hancox et al., 2014).Among debris (colluvium) failures, the failure mode will vary, depending on the material content and whether failure took place in brittle or ductile zones.

Production of landslide inventories
In order to produce regional inventories of landslides triggered by these events, landslide scars were identified and mapped through stereoscopic interpretation of panchromatic aerial photographs, combined with ground and oblique aerial photography, based on morphometric criteria and the surface reflectivity contrasts between undisturbed ground and failures (Nichol and Wong, 2005;Liu et al., 2002;Hovius et al., 1997).Our inventories consist of landslides where the failed mass evacuated the failure plane and moved downslope to leave a discernable, bare-earth scar.Accordingly, all landslides included in the inventories are disrupted slides, which moved rapidly downslope following failure.
Landslides triggered by the 1929 earthquake were mapped using 1 : 86 000 scale images taken in February 1968, and validated using ground photos taken in 1929 and further aerial photos taken in 1947 (SN 265, runs 1457(SN 265, runs -1463) ) for selected regions (Appendix A).From comparison of earlier and later imagery, we found that scars from landslides triggered in 1929 were still clearly visible and could be mapped in imagery acquired 39 years after the earthquake, due to a slow rate of regeneration of native bush.This is particularly true for larger, bedrock failures, while smaller soil and de- bris failures are more rapidly obscured by vegetation.Landslides attributed to the 1968 earthquake were mapped using 1 : 66 000 scale panchromatic aerial images (Appendix A) taken in November 1974 and aerial oblique and ground photos taken in 1968-1969.Landslide mapping was further validated based on observations from fieldwork undertaken by G. Hancox throughout 1968 and1969, and during aerial reconnaissance undertaken by G. Hancox in 1998and 2010, and in 2011 by R. Parker.Comparison of pre-and post-1968 imagery was carried out to delineate 1929 landslide areas from those triggered or further influenced by the 1968 earthquake (Fig. 7).Although the intervening periods between seismic events and imagery acquisition create potential for inclusion of landslides triggered by aseismic (rainfall) events, observations from reconnaissance between 1968 and 2014 and historical records compiled by the West Coast Regional Council (Hancox et al., 2014) suggest a lack of widespread landsliding resulting from heavy rainstorms or other processes during inter-seismic periods, supporting a seismic mode of triggering for the landslides observed (Pearce and Oloughlin, 1985;Hancox et al., 2014).Prior to the 1929 earthquake, two large events of M w ∼ 7 are estimated to have occurred in 1868 and 1893, with epicentres located around 200 km to the north and northeast of the study area, respectively (Anderson et al., 1994).Due to the lack of pre-1929 imagery, there may be potential for landslides triggered by these events to be wrongly attributed to the 1929 earthquake.However, due to their distance from the study area, these events would have produced relatively weaker ground motions than the 1929 event -MMI V-VII (1868, 1893) vs. MMI IX-X (1929) (Anderson et al., 1994;Hancox et al., 2002) -capable of triggering few, relatively small landslides (Hancox et al., 2002).As smaller landslides are more rapidly obscured by vegetation, it is unlikely that smaller failures from these events feature in our data set.An earlier larger earthquake of around M w ∼ 7.4 is also thought to have occurred ca.1650, as indicated by  several landslide-dammed lakes in the northwest Nelson area (Hancox et al., 2002;Perrin and Hancox, 1992;Henderson, 1937).Larger, visible pre-20th century landslide scars in the region were mapped separately and are not included in this analysis.
Polygons delineating the combined landslide source and runout areas of individual landslides were mapped by hand on 1 : 50 000 scale topographic maps, which were then digitised and imported into a GIS.Particular effort was made to map individual failures separately and separate coalesced landslide features, in order to avoid issues of feature amalgamation in the data set (Li et al., 2014).The imagery resolution allowed mapping of landslides down to a minimum size of ∼ 50 × 50 m (∼ 2500 m 2 ).For the 1929 earthquake, 4074 landslides (182 km 2 total landslide area) were mapped across an area of 4222 km 2 .Note that this mapping covers the southern half of the landslide-affected area, while the 1929 landslides extend to the north, away from the region affected by the 1968 earthquake.By contrast, for the 1968 Inangahua earthquake, 1400 landslides (39 km 2 total landslide area) were mapped across an area of ∼ 3500 km 2 .Of these, 246 landslides were reactivations or enlargements of landslide scars that failed in 1929, mostly in over-steepened source areas of the previous failures.The areal extents of the landslide inventories overlap by 2882 km 2 , ∼ 80 % of which experienced MMI ≥ VIII in both events.The areas of both the 1929 and 1968 landslides exhibit characteristic powerlaw scaling (e.g.Hovius et al., 1997;Guzzetti et al., 2002;Malamud et al., 2004;Van Den Eeckhaut et al., 2007;Fig. 8): where p(x) is the probability of a landslide having a given area, x min is the minimum area of landslide modelled by the  Meunier et al. (2007Meunier et al. ( , 2013)), Dai et al. (2011), Lee et al. (2008), Hancox et al. (1997Hancox et al. ( , 2002) ) Arias intensity, MMI) Distance from the seismic source Regional attenuation of seismic wave amplitudes Position on hillslope (normalised Ridge to stream patterns of Davis and West (1973), Bouchon (1973), Wu et al. (1990), Benites et al. (1994), Meunier et al. (2008), Densmore et al. (1997) distance from stream to ridge crest) topographic amplification and damping Orientation of hillslope relative Directional patterns of topographic to seismic source amplification and damping, due to the incidence angle of seismic waves Hanging wall vs. footwall location Proximity of the fault and enhanced Abrahamson and Somerville (1996), of sites rupture directivity effects in hanging wall areas Somerville et al. (1997), Abrahamson et al. (2008) Strength of hillslope materials

Rainfall
The effect of pore water pressure in reducing Dellow and Hancox (2006), hillslope effective stress Iverson (2000) Static stress loading in hillslopes Hilllope gradient Magnitude of static stress loading in hillslopes Keefer (2000), Khazai and Sitar (2004), Lee et al. (2008), Dai et al. (2011) Local hillslope relief function and α is the power-law scaling exponent.The positions of the rollover for smaller landslides suggest complete mapping of landslides larger than between 11 000 and 13 000 m 2 in both data sets.These limits result from the combined effects of the mapping resolution and vegetation recovery in obscuring smaller landslide scars.The power-law scaling exponents of 2.68 (1929) and 2.85 (1968), fitted using the method of Clauset et al. (2009), fall within the typical range of previously observed values for landslide inventories (1.4-3.4), which have a central tendency around 2.3-2.5 (Van Den Eeckhaut et al., 2007;Stark and Guzzetti, 2009).
To analyse the spatial pattern of hillslope failures, we use the landslide source areas, rather than areas covered by landslide runout and deposits.For most landslides it was difficult to visually separate landslide source and runout or deposit area.Based on a sample of 51 landslides where visual delineation of the source area was possible, dividing the extent of each landslide at its midpoint elevation (i.e. the con-tour halfway been the maximum and minimum landslide elevation) provided a good approximation of the separation between source and runout-deposit (Appendix B).This approach is similar to the method of extracting landslide areas above the median landslide elevation, which has been employed in previous studies (Parise and Jibson, 2000;Jibson et al., 2000;Capolongo et al., 2002;Lee et al., 2012).However, our technique is less prone to overestimation of the source area for landslide masses that run out over large distances across low-gradient ground.The residuals are calculated by aggregating probabilities across equal quantile bins of the x variable.Note that positive and negative residuals are relative to the prediction of a model that does not explicitly consider the effect of hillslope preconditioning but is fitted using landslide data that are subject to the effect of hillslope preconditioning.Therefore it is the direction of the trend, rather than the absolute (positive or negative) residual values, that is of importance in the test of preconditioning.
that influence the intensity of event-specific seismic ground motions, and the strength of hillslope materials and the static shear stresses acting on them, which may remain more consistent from one event to the next.Empirical studies have revealed a number of proxy variables that can be used to represent these factors at the regional scale (Table 2).
Logistic regression is a standard technique for assessing controls on earthquake-triggered landslide distributions (e.g.Yesilnacar and Topal, 2005;Dai and Lee, 2003;Garcia-Rodriguez et al., 2008;von Ruette et al., 2011), by modelling the influence of multiple predictor variables on a categorical response (Cox, 1958;Walker and Duncan, 1967).The function takes the form where logistic regression is used to estimate the coefficients (b, b n . ..) for predicting the probability that Y = 1, given the values of one or more predictor variables (x, x n . ..).In this case, Y = 1 corresponds to the occurrence of a landslide at a particular point in space.
Although previous studies have applied logistic regression with the implicit assumption of temporally static hillslope sensitivity to landslide triggering, here we use this technique to test a hypothesis of hillslope preconditioning for failure by previous events.We first undertake an implicitly static logistic regression analysis in order to model the distributions of landslides, as can best be achieved without considering the influence of past events.We hypothesise that if the 1929 earthquake influences the 1968 landslide distribution, then the residual variability, unexplained by our regression model, must exhibit a relationship with the spatial distribution of the effect of the previous earthquake on hillslopes.To test this hypothesis, we compare the residuals of our 1968 regression with a measure of hillslope preconditioning, here the probability of landslide occurrence in 1929.A graphical representation of hypothetical outcomes is presented in Fig. 9.We assume that logistic regression models have been fitted and used to hindcast the probability of hillslope failure (P LS ) for both earthquakes.Note that by definition the observed probability of landsliding, being based on observations a posteriori, is 1 for landslide sites and 0 for non-landslide sites.For comparison, observed and predicted probabilities are therefore aggregated (mean-averaged) across sites (pixels) that fall within equal quantile bins of the predictor variables.For each data point generated, the mean predicted probability represents the proportion of sites expected to fail, while the mean observed probability represents the proportion of sites observed to fail.If the model for the 1968 earthquake is accurate, then the residuals (observed P LS minus predicted P LS ) should yield no structure when plotted against the predicted values (Fig. 9a).Similarly, there should be no structure in the residuals when plotted against each of the individual predictor variables (Fig. 9b).However, if the 1929 earthquake has influenced the 1968 landslide distribution, then the residuals should exhibit structure when plotted against the predicted P LS for the 1929 earthquake.Fig. 9c illustrates two endmember scenarios, showing how the 1929 earthquake might be expected to influence the 1968 landslide distribution: 1. Hillslopes with higher predicted P LS in 1929 exhibit lower than expected P LS in 1968.This could be the case if widespread failure of unstable hillslopes in 1929 resulted in fewer hillslopes being "available" for failure in 1968.
2. Hillslopes with higher predicted P LS in 1929 exhibit higher than expected P LS in 1968.This could be the case if, despite the widespread failure of unstable hillslopes in 1929, damage accumulation in those hillslopes that did not fail primed those sites for failure in 1968.
Conversely, if there were no trend in the residuals, this would suggest that the 1929 earthquake has not influenced the 1968 landslide distribution.Although damage accumulation is specific to landslides in brittle hillslope materials, and not necessarily present in all hillslopes where landslides have been mapped, even if a subset of hillslopes record the legacy of past earthquakes, we should expect to see the signal via this test.
In order to undertake logistic regression analysis, we first removed landslides with areas less than 13 000 m 2 from our data set to eliminate biases arising from small landslides censored by the mapping resolution and post-landslide vegetation regrowth.We then defined a sample grid at 30 m resolution, based upon a digital elevation model, resampled from the 10 m resolution New Zealand Digital Terrain Model (GNS Science, 2011), using bilinear resampling.The elevation model was resampled at this scale to remove fine-scale noise, while ensuring that the characteristics of individual landslides are resolved.Using a 30 m grid, we ensure that more than 10 sample points fall within the smallest landslides included in our analysis.Additionally, 30 m is much less than typical hillslope lengths in the region of 500 m, ensuring that multiple hillslopes are not contained in a single pixel.Response and predictor variables were then generated for each grid cell.For the response variable, binary grids of landslidesource and non-landslide-source pixels were generated from the mapped 1929 and 1968 landslide source zones.We removed from this analysis the 246 landslides from the 1968 data set that occurred as reactivations of 1929 landslide scars.This allowed our analysis to test exclusively for the influence of hillslope damage accumulation, rather than the effect of slopes over-steepened or undermined by previous landslides.Predictor variables (Fig. 10, Table 3) were derived to represent factors previously found to influence landslide occurrence elsewhere (Table 2).For both earthquakes, we used the horizontal distance of each grid cell to the surface projection of the fault (FLD), and the three-dimensional distance from each grid cell to the closest point on the coseismic fault plane (FPD) as proxies for the regional attenuation of seismic waves and shaking intensity.For the 1968 earthquake, we also used the Shakemap PGA model for this purpose by interpolating from modelled PGA values at 0.05 • (∼ 4.5 km) grid spacing (PGA).A binary variable, HW, coding the hanging walls (HW = 1) and footwalls (HW = 0), was used to Each map shows distributed values of each predictor variable across the 5629 km 2 combined area of landslide mapping for both events (as shown in Fig. 1).Variable descriptors and units are summarised in Table 3.
represent hanging wall effects on ground motion (Abrahamson and Somerville, 1996).A second binary variable, DIR, coding regions towards and away from which the fault ruptures propagated (0 and 1, respectively), was used to represent the effect of rupture directivity on ground motion (Bray and Rodriguez-Marek, 2004).The local hillslope orientation (HO) relative to the seismic source (0 for hillslopes with aspect oriented away from the fault rupture, and 1 for aspects oriented towards the fault rupture) was used to represent the incidence angle of seismic waves.Normalised distance from stream to ridge crest (0 for sites located in a stream channel, 1 for sites located on a ridge crest) was used to represent valley-scale patterns of topographic amplification and damping (NDS).Local hillslope gradient, measured over a three-pixel (90 m) spatial window (SL) and two relief metrics (the relief (ER) and standard deviation (ES) of elevation within individual drainage basins, divided by the drainage Binary variable coding the 1929 hanging wall and footwall HO-1929 Local hillslope orientation relative to the 1929 seismic source (incidence angle of seismic waves) •

FLD-1968
Horizontal distance of each sample cell to the surface projection of the 1968 fault km FPD-1968 3-dimensional distance from each sample cell to the closest point on the 1968 coseismic fault plane km HW-1968 Binary variable coding the 1968 hanging wall and footwall HO-1968 Local hillslope orientation relative to the 1968 seismic source (incidence angle of seismic waves) Long-term mean antecedent precipitation total for 3 months prior to the earthquake mm PD6 Long-term mean antecedent precipitation total for 6 months prior to the earthquake mm basin area) were used to represent the magnitude of static stresses.In the calculation of elevation derivatives, using a spatial window size (three pixels or 90 × 90 m) smaller than the smallest individual landslides included in our analysis, we minimise the risk of overgeneralising the characteristics of individual landslides.A categorical variable indicating different lithologies was used to represent variability in material strength (G).In order to capture the regional distribution of structure on bedrock landslides, we generated a binary variable of dip/anti-dip slopes (DS) by comparing local slope gradient and aspect with the azimuth and dip of recorded structures from the New Zealand QMap data set (Rattenbury et al., 1998(Rattenbury et al., , 2006;;Nathan et al., 2002), which was interpolated using Thiessen polygons.The northerly component of aspect (cosine of aspect, CA) is used to characterise hillslope-scale variations in received solar radiation, which have been associated with the relative intensity of physical and chemical weathering (Mcfadden et al., 2005).Note that CA = 1 indicates north-facing hillslopes, which experience higher levels of southern hemisphere solar radiation, while CA = −1 indicates south-facing hillslopes.
Although a clear relationship between rainfall-induced pore pressure and earthquake-induced landslide triggering has not be shown, collated data sets for New Zealand suggest that earthquakes that occur during wetter months trigger more landslides than those during drier periods (Dellow and Hancox, 2006).Rainfall records from Karamea (see http://cliflo.niwa.co.nz/) suggest similar levels of rainfall preceded the two events, which had a difference in calendar position of only 23 days.For example, June 1929 received 307 mm (May-June 1929 received 406 mm) and May 1968 received 275 mm (April-May 1968 received 525 mm).It is therefore unlikely that a difference in groundwater conditions is responsibly for a difference in sizes of landslide events triggered by the earthquakes.On the scale of the individual events, we account for spatial variability in rainfall across the region using mean monthly precipitation totals for the period 1950-2000 (Hijmans et al., 2005).These were used to estimate antecedent precipitation totals for each grid cell, for the 3-month (PD3) and 6-month (PD6) periods prior to each earthquake.Note that these data are used to indicate the generalised mean distribution of rainfall across the region, rather than data for those specific years.In the absence of high-resolution spatial-temporal rainfall data, this allows our analysis to test whether regions of the landscape receiving more rainfall are more susceptible to landslide triggering than drier regions.
In order to avoid the problem of over-fitting regression models and predictor covariance, issues particularly characteristic of automated fitting procedures (e.g.Hosmer and Lemeshow, 2000), model fitting was undertaken manually, and based on the following criteria: 1.All predictors must have a logical, statistically significant (p < 0.05) and consistent influence on P LS for both earthquakes.Whilst the regression coefficient associated with a variable may differ between the two events, this condition stipulates that the direction of influence (+/−) must remain constant.
2. Predictor variables included in the model must not exhibit multicollinearity, as determined by variance inflation factors (VIF): where ln L(M full ) is the log likelihood of the full model and ln L(M intercept ) is the log likelihood of the model without any predictors.The pseudo-R 2 is designed to look like a conventional R 2 goodness of fit, derived from ordinary least square regression, with values ranging from 0 (no correlation) to 1 (perfect correlation or, in the case of logistic regression, perfect separation of true (landslide) and false (non-landslide) categories).As logistic regression is fitted through an iterative process of maximum-likelihood estimates, the conventional R 2 approach to goodness of fit does not apply.However, like conventional R 2 values, pseudo-R 2 can be seen as an indicator of explained variability and the level of improvement offered by the full model over the model without its predictors (McFadden, 1974).
During the fitting process, multiple variable combinations were iteratively tested.The final models presented below represent those that produced the best fit whilst meeting the above criteria.
During the model fitting, grid cells with hillslope gradient > 58 • were found to produce numerical problems associated with the very low frequency of data at high values.This amounted to an area of 1.3 km 2 (less than 0.05 % of the study area).In this range the relationship between hillslope gradient and failure probability was found to exhibit a rollover, suggesting a decrease in failure probability at high gradients.It is unclear whether this behaviour is real, an artefact of the low data frequency, a reflection of the difficulty of mapping landslides on steep slopes from aerial imagery, or a deterioration of DEM quality at high gradients.As the logistic function cannot simulate a modal (humped) relationship, and as slope gradient is one of the dominant variables in the model, these cells were removed from the analysis prior to model fitting.5) (fault distance model) and (e) output map of predicted P LS from Eq. ( 6) (PGA model).Plots of observed vs. predicted P LS : (f) 1929 earthquake, Eq. ( 5); (g) 1968 earthquake, Eq. ( 5); and (h) 1968 earthquake, Eq. ( 6).These data are generated by aggregating probabilities across 20 equal quantile bins of the predicted P LS .

Earthquake-induced hillslope failure probability models
We derived two fitted model versions to hindcast hillslope failure probability; these differ in their characterisation of the regional distribution of ground motions.For both earthquakes, models were derived using fault plane distance, FPD, and location relative to rupture directivity, DIR, as a proxies for ground motion.For the 1968 earthquake we also present a model using PGA in place of FPD + DIR, which constrains the landslide distribution more accurately.In our FPD-based model for the 1929 and 1968 earthquakes, hillslope failure probability can be modelled via the following equation: where the regression coefficients are indicated by c.Similarly, in our PGA-based model for the 1968 earthquake, hillslope failure probability can be modelled via the following equation: The regression coefficients and fit statistics for these models are given in Table 4, while Fig. 11 presents a comparison of predicted and observed P LS .
In each model, landslide probability is expressed as a function of the regional seismic ground motion (characterised by three-dimensional distance from the fault plane and location relative to rupture directivity, or Shakemap PGA), hillslope gradient (where the influence of hillslope gradient varies with lithology) and normalised distance from stream to ridge crest.Note that the lithology predictor variable has more explanatory power and significance when it is used to allow variability in the effect (coefficient) for hillslope gradient, rather than allowing a categorical lithology variable to modify landslide probability directly.All other variables tested during model fitting were found to be less effective predictors than those included in the models presented, or they failed in one or both of the fitting criteria.Note that these models describe the relative spatial distribution of landslides, while absolute differences in the magnitude of the earthquakes are accounted for implicitly by fitting the model separately for each earthquake.
For both model versions, predicted and observed probabilities display a good fit to the line of equality (Fig. 11f-h).P LS values hindcast using Eq. ( 5) display a slight over-prediction at low probability values.It is likely that these errors at the lower limit of the distribution of probabilities are at least in part statistical artefacts of low data frequency and near-zero probability values.Spatially, values hindcast using Eq. ( 5) also display a step change produced by the DIR binary variable (Fig. 11b and d).Although this discrete artefact is unrealistic, the variable exhibits a significant, physically plausible effect, which improves the fit of both models; landslide probabilities are higher in regions towards which the rupture propagated, where stronger ground motions are expected (Bray and Rodriguez-Marek, 2004).P LS values hindcast using Eq. ( 6) do not exhibit this artefact (Fig. 11e), and provide a better statistical fit to the observed landslide distribution.
Although each predictor explains a component of the variance in the spatial distribution of failures, not all variables  5) for the 1929 earthquake, (b) results from Eq. ( 5) for the 1968 earthquake, and (c) results from Eq. ( 6) for the 1968 earthquake.contribute equally.Figure 12 presents the predictor variables in rank order of their importance in each model, determined by sequentially removing the predictor contributing least to the fit of the model.In all three models, the regional ground motion proxy (distance from the fault plane or PGA) and hillslope gradient rank as the top two variables, followed by geology, position on hillslope, and location relative to rupture directivity in the case of Eq. ( 5).In all three models, the regional ground motion, hillslope gradient and geology account for over 80 % of the total model fit, while position on hillslope (NDS) and location relative to rupture directivity (DIR) are secondary in defining the spatial distribution of landslides.The position on hillslope (NDS) relationship is consistent with ridge-to-valley scale patterns of amplification and damping of seismic waves found by others (Davis and West, 1973;Bouchon, 1973;Wu et al., 1990;Benites et al., 1994;Meunier et al., 2008).
Note that in none of the models does predicted P LS ever reach values of 1 or 0, as the predictor variables are not able to discriminate slopes where failure or non-failure is a certainty.This observation may be attributed to stochastic uncertainty in the model predictors, but it also points to the possibility of important factors omitted from the model (epistemic uncertainty), of which the unconstrained damage legacy of past events is one possible candidate.While we do not have data to constrain the influence of all past events that have possibly conditioned hillslope materials, we are now able to test whether the damage legacy of the largest recent earthquake (the 1929 earthquake) may be present in the distribution of landsliding triggered by the 1968 earthquake.Positive residuals indicate that the model under-predicts P LS and negative residuals indicate that the model over-predicts P LS .Note that the amplitude of the plotted residuals varies due to binning and aggregating probabilities with the different predictor variables.For each plot, the coefficient (r) and significance (p) of correlation in the residuals are given.These were derived by adding each x variable into a logistic regression analysis of P LS predicted using Eq. ( 6) and observed landsliding.In this respect we test the coupled significance of SL and G as they feature in the model.

Potential influence of the 1929 earthquake on the 1968 landslide distribution
We use Eq. ( 6) to test our hypothesis of the influence of the 1929 earthquake on the landsliding resulting from the 1968 event.This model most accurately hindcasts the 1968 landslide distribution, using a ground motion term -PGA (1968) -that cannot be overfitted to the landslide distribution.Conversely, in Eq. ( 5), DIR forms part of the ground motion term, giving the model an additional degree of freedom to account for any imbalance in P LS between the northeast and southwest quadrants.As the northeast quadrant represents much of the area closest to the 1929 source, while the southeast quadrant is further away, the DIR variable absorbs and masks some of the preconditioning signal of the 1929 earthquake.
Any test for preconditioning depends on having a regional ground motion term that cannot be overfitted in this way, for which we use the Shakemap PGA field.As this is based on observed ground motions that are influenced by rupture directivity, the Shakemap PGA also implicitly accounts for the effect represented by DIR in Eq. ( 5).However, as the Shakemap data are subject to large uncertainties, we present the following result as tentative, using the best available data for these events.To test whether the 1929 earthquake has influenced the 1968 landslide distribution, we use FPD (1929) as a proxy of the regional distribution of ground motion pro-duced by the 1929 earthquake, in the absence of PGA data for this event.We acknowledge that this scenario is not ideal and it would be preferable if PGA (1929) were available, along with PGA (1968).However, in the absence of 1929 PGA data, a distance term provides a reasonable proxy for the spatial pattern of ground motions (Campbell and Bozorgnia, 2008).Figure 13 presents the results of our analysis in a form equivalent to that outlined conceptually in Fig. 9, with correlation coefficients (r) and p values to test the strength and significance of trends in the residuals.When tested against the predictor variables (Fig. 13a-d), there is no monotonic trend and little structure in the residuals, which suggests that the model predictors are well fitted to the data.There is minor nonlinearity in the residuals plotted against NDS, which suggests that the increase in P LS with NDS begins to saturate close to the top of hillslopes.While this results in slight over-prediction of near-ridge-top landslide probability (sites with NDS > 0.7), we found that removing these sites from our analysis does not change the result that we now discuss.
When plotted by hindcast P LS for the 1929 earthquake (P LS 1929), the 1968 residuals display a significant positive trend (Fig. 13g).Hillslopes with P LS 1929 greater than 0.013 (38 % of the overlap region mapped for both events) exhibit higher P LS in the 1968 earthquake than predicted Earth Surf.Dynam., 3, 501-525, 2015 www.earth-surf-dynam.net/3/501/2015/FPD (1929) into Eq.( 6), we are able to improve the fit of the logistic regression model from R 2 = 0.246 to R 2 = 0.251 (Table 5).Additionally, by adding DIR (1968) together with FPD (1929) into Eq.( 6), we can check whether this result can be attributed to a lack of consideration given to rupture directivity in the Shakemap PGA data (Table 5).Although the rupture directivity term has a significant coefficient, the relationship between P LS and FPD (1929) is still present and statistically significant.

Discussion
Our analysis has sought to control for all major factors known to influence the spatial distribution of landslides, at (or close to) the scale of the whole earthquake-induced landslide event.Using the best available data for the 1968 earthquake, our model achieves this using variables with defined physical links to landsliding while maintaining a low level of model complexity, which avoids overfitting.Once these steps have been taken to control for the influence of other variables, our results suggest that landslide probability in 1968 is higher for hillslopes that experienced strong ground motions in the previous 1929 earthquake.
Our results support the findings of previous work into modelling earthquake-induced landslides, as well as providing new insights into how past earthquakes may influence future landslide distributions.The roles of individual components in our logistic regression models are in agreement with those observed in previous studies (e.g.Dai et al., 2011;Meunier et al., 2008;Meunier et al., 2007), and the presence of these relationships in both the 1929 and 1968 earthquakes supports the extrapolation of these models both temporally and spatially.A number of variables that we might expect to influence the landslide distribution showed no significant influence when the effect of other predictors was controlled for.This particularly concerns factors influencing the aspect of landslides.Neither the orientation of hillslopes relative to the seismic source nor relative to hillslope-scale variations in received solar radiation was found to exhibit a significant influence on landslide probability.This implies that patterns observed in other earthquakes may be regionally specific or confounded by the influence of other more "powerful" predictors that might not have been controlled for and which co-vary with the aspect of hillslopes.For example, given the dominant influence of hillslope gradient on landsliding, systematic variations in hillslope gradient with aspect are a possible candidate for this confounding effect.
The consistency with which the model describes the spatial distribution of hillslope failures for both events suggests that the combination of underlying relationships presented in Eqs. ( 5) and ( 6) may be applied more generally to earthquakes in this region.In other words, landslides triggered by earthquakes in this area are likely to conform to the spatial distribution of hillslope failure probability described here.By removing the less influential variables and identifying the major regional-scale influences on failure probability, the model can be made less event-specific and therefore more transferrable.The combination of ground motion and local hillslope gradient, with the influence of hillslope gradient dependent on lithology, therefore provides a candidate variable subset for a generalised earthquake-induced landslide probability model.
While time-independent variables provide useful constraints on the spatial distribution of earthquake-triggered landslides, our results also suggest that previous earthquakes may impart an influence on future landsliding.Residuals in the landslide distribution predicted for the 1968 earthquake suggest that hillslopes with higher predicted P LS in 1929 (or those closer to the 1929 seismic source) exhibit higher than expected P LS in 1968.This implies that, despite the widespread failure of unstable hillslopes in 1929, at least some of those hillslopes that did not fail in or shortly following 1929 were more susceptible to failure in 1968.This behaviour is consistent with our hypothesised influence of damage accumulation, where failure occurred in brittle hillslope materials.Our results suggest the possibility that in the case of the 1929 earthquake, damage in unfailed hillslopes persists, resulting in regions close to the 1929 seismic source having enhanced sensitivity to landslide triggering in 1968.We stress that this suggestion must be treated as tentative due to uncertainties in our analysis variables.This particularly applies to the ground motion proxies and the PGA field, which relies on interpolation from limited observations, using ground motion prediction equations (Wald et al., 2006).Additionally, as the elevation model used in our analysis was derived following both earthquakes, there is the possibility that hillslope gradients measured at landslide sites may not accurately reflect slope characteristics at the time of landslide triggering.However, depths of most mapped landslides are likely to be smaller than uncertainties in the elevation data, suggesting that the 1929 and 1968 landslides are unlikely to have produced surface changes detectable in the elevation model (Appendix D).As our analysis explicitly considered only the source area of landslides, any bias is likely to involve overestimation of gradients, in source areas where head scarps have been steepened by landsliding, or no effect in cases of translational failures.In our probability modelling, underestimation of gradients for landslide sites produces over-prediction of landslide probability for steeper hillslopes.The residuals of Eq. ( 6), plotted against slope gradient (Fig. 13b), may indicate very slight over-prediction at high gradients, reflecting this effect.However, as postlandslide topographic changes are small relative to the elevation model uncertainty, and as slope gradient appears to be well fitted to the data, this suggests that the use of postlandslide elevation data is unlikely to effect the outcome of our analysis.
Further advances in testing our theory may be made where multi-earthquake landslide data sets are available for more recent events, where higher resolution (and multi-temporal) elevation models are available, along with data from more dense seismic networks.On the basis that future testing may further support our hypothesis, we discuss the implications of our results in light of current understanding of the temporal landslide response to earthquakes.Fundamentally, our results are consistent with the idea that seismic ground motion produces irreversible damage, such that the legacy of past earthquakes may be preserved to a greater or lesser degree as a loss in strength in hillslope materials, for longer periods of time than previously thought.Several studies suggests that, following large earthquakes, prolonged rates of mass wasting, and associated indicators of changes in hillslope material strength, return to background levels within timescales of less than a decade (Hovius et al., 2011;Uchida et al., 2014;Marc et al., 2014).However, our data suggest that even after several decades, when the next large earthquake occurs, there is still a signal of hillslopes weakened by the previous earthquake.Note that, unless some healing or annealing process takes place in hillslope materials, or all damaged material is stripped from hillslopes by erosion, there is no reason why we should not expect this to be the case.We explain this observation further by considering groups of hillslopes in different states from a spectrum of earthquake-induced damage.During the 1929 earthquake a first subset of hillslopes is weakened to the point of failure (co-seismic landslides).A second subset of hillslopes is moved to states close to the point of failure, such that failure of these hillslopes is triggered during relatively moderate aftershocks and post-seismic rainfall events.Landslides produced by these two subsets of hillslopes generate sediments that take time to be evacuated from the orogen by fluvial processes, at a rate that decays over a sub-decadal timescale as landslide deposits are exhausted of mobilisable sediment (Hovius et al., 2011;Dadson et al., 2004).A third subset of hillslopes has also been weakened by the 1929 earthquake, but insufficiently for moderate post-seismic events (aftershocks and rainstorms) to trigger failure.Additionally, yield stresses in these hillslopes may remain too high to be exceeded by moderate interseismic events, such that continued permanent deformation and damage accumulation does not occur post-seismically.As a result, the post-seismic rate of landsliding decays, while the landscape maintains a subset of hillslopes damaged and in a state closer to failure than prior to the 1929 earthquake, but which may only be brought to the point of failure by another large earthquake.Both coseismically and post-seismically, only a relatively small proportion of hillslopes in the landscape actually undergo full failure.For example, within 10 km of the 1929 source, only 3 % of hillslopes were mapped as 1929 landslides.Therefore the behaviour of hillslopes that fail during or soon after an earthquake only accounts for small subset of the landscape effected by the seismic ground motions.The result we present here, and numerical simulations using geotechnical models (Parker et al., 2013;Moore et al., 2012), supports the hypothesis that there is a legacy of damage in the remaining apparently intact landscape that may not fail either during or after an earthquake.If this is the case, then, at any point in time, each of these subsets exists along a continuum from pristine hillslopes to those damaged almost to the point of failure, evolving with each event that generates damageinducing stresses.
This long-term perspective may reveal why correlation between the 1968 landslide and the 1929 earthquake is weak.Although our analysis provides spatial estimates of the effect of the 1929 and 1968 earthquakes on hillslopes, we lack information on the damage condition of hillslopes prior to the 1929 earthquake.Hence we can only expect to find partial or weak correlation with a single past event, even if the 1968 landslide distribution were the deterministic product of the accumulation of all past events.However, one would expect events added to the historical record to incrementally and cumulatively account for more unexplained variability in landsliding.
Similarly, if landslide distributions are predetermined by the legacy of accumulated damage from past events, then data from neither the triggering earthquake nor a single previous event can provide an exact prediction of landsliding.In this way, the apparently stochastic nature of landslide occurrence and the inability of current models to identify the exact hillslopes that undergo failure may in part result from not knowing the condition of each hillslope at the onset of shaking.In future, if the damage condition of hillslopes can be correlated with the history of past damage-inducing events, then building historical data or proxies for damage into landslide models may provide a means of constraining this effect.

Conclusions
The main conclusions drawn from this study can be summarised as follows: 1.The 1929 and 1968 earthquakes reveal a consistent spatial pattern of landslides that can be modelled probabilistically as a function of spatial variability in seismic ground motion, hillslope gradient, lithology and position on hillslope.Statistically, the seismic ground motion and hillslope gradient (where the influence of hillslope gradient is lithologically dependent) account for the majority (> 80 %) of the explanatory power of the model.We therefore conclude that these factors are the most important considerations for predicting an earthquake-induced landslide distribution at the regional or whole-event scale.
2. Our results suggest that the legacy of the 1929 Buller earthquake may have influenced the spatial distribution of landslides triggered by the 1968 Inangahua earthquake, in a manner consistent with the accumulation of damage in hillslopes that did not fail in the 1929 earthquake.Although uncertainty in input variables makes our result necessarily tentative, our methodology could be used for further testing of this hypothesis where multi-earthquake landslide data exist for more recent events.
3. The signal of the 1929 earthquake revealed in landslides triggered by an earthquake several decades later suggests that the damage legacy of past earthquakes persists in the landscape for much longer than suggested subdecadal decay periods of post-seismic landslide activity.We speculate that this may be due to damage that persists in hillslopes that do not fail co-or post-seismically but have only been sufficiently weakened to fail in the next large triggering event.
4. Long-term damage accumulation provides a potential deterministic explanation for the observed stochastic spatial nature of landslide triggering.Predetermined by the legacy of damage accumulated in hillslopes during past events, landslide distributions may be viewed as a function of the damage history of the landscape surface.
We may expect each past event to partially correlate with unexplained variability in landsliding, resulting in improved predictions as the historical record of triggering events is extended.
5. Accurate, multi-event mapping of landslide distributions, resampled following individual earthquakes and storms, may represent a significant step towards better understanding of temporal correlation between past and future landslide-triggering events.To implement the findings of this investigation into regional landslide assessments, future work should explore the value of adding historical and palaeo-seismic and climatic data into landslide models, providing a means of making susceptibility assessments dynamic through time.

Figure 1 .
Figure 1.Elevation map of the Buller River to Karamea of the northwest South Island, New Zealand, showing the sources, ground motions and landslides triggered by the 1929 Buller and 1968 Inangahua earthquakes.Included on this map are the epicentres, focal mechanisms and fault planes of the earthquakes (Ander-son et al., 1993, 1994;Stirling et al., 2007), isoseismal contours(Dowrick, 1994;Adams et al., 1968), mapping coverage regions and earthquake-induced landslides mapped in this study.

Figure 2 .
Figure 2. (a) 1968 rock falls from limestone bluffs at White Cliffs Escarpment, 5 km west of Inangahua.(b) 1968 rock and debris fall area in the upper Buller Gorge.

Figure 3 .
Figure 3. (a) 1968 debris slide on road cut in the upper Buller Gorge.(b) Multiple 1968 debris slides on slopes on the south side of the Buller River.

Figure 4 .
Figure 4. (a) Lake Stanley rock avalanche.The lake was dammed by the landslide, which was triggered by the 1929 earthquake.(b) 1995 photo of a rock avalanche that dammed the Buller River during the 1968 earthquake.Apart from vegetation growth, the scar has changed little in the last 40 years.

Figure 5 .
Figure 5. (a) Matakitaki landslide triggered by the 1929 earthquake.Debris from this large (18 million m 3 ) dip-slope rockslide travelled ∼ 1 km across the valley floor, killing four people and forming a landslide dam.Note that, after over 70 years, the landslide scar is still visible.(b) Aerial view of the Matiri Valley (15 km north of Murchison), which was extensively damaged by landslides during the 1929 earthquake.Numerous scars of rockfall and debris slides are still clearly visible in 2011.

Figure 6 .
Figure 6.(a) 1968 Oweka rock slide with rock debris (Lensen and Suggate, 1968).(b) 1968 rotational slide of ∼ 2 million m 3 in sandy ("Blue Bottom") mudstone.At the top of the landslide the semiintact block below the prominent head scarp has slumped about 6 m.The main body of the slide has carried the rock downslope and comprises highly disrupted mudstone boulders and finer debris.

Figure 7 .
Figure 7. Two aerial images used to map landslides triggered by the 1929 and 1968 earthquakes in the Buller Gorge area.Scars from the 1929 landslides are recognisable on the Survey No. 2033 image (a), 39 years after the earthquake, and many scars were reactivated or enlarged by the 1968 earthquake (b).

Figure 8 .
Figure 8. Probability densities of landslide areas as a function of landslide area for the 1929 Buller and 1968 Inangahua earthquake landslide inventories.Power-law scaling exponents(α) have been derived using the method of Clauset et al. (2009), for areas > 11 000 m 2 .

Figure 9 .
Figure 9. Hypothetical output of hillslope failure conditional probability analysis for the 1929 and 1968 earthquakes.This illustrates how the model residuals (observed P LS minus predicted P LS ) would be distributed if conditional probability models, which account for the full subset of spatial, time-independent factors (i.e.those not associated with previous events) influencing landslide occurrence, have been fitted and used to hindcast P LS .(a) 1968 model residuals are uncorrelated with predicted P LS .(b) Model residuals are uncorrelated with values of individual model predictors.(c) Model residuals are correlated with P LS hindcast for the 1929 earthquake either negatively (model 1)indicating preconditioning of hillslopes against failure -or positively (model 2) -indicating preconditioning of hillslopes for failure.The residuals are calculated by aggregating probabilities across equal quantile bins of the x variable.Note that positive and negative residuals are relative to the prediction of a model that does not explicitly consider the effect of hillslope preconditioning but is fitted using landslide data that are subject to the effect of hillslope preconditioning.Therefore it is the direction of the trend, rather than the absolute (positive or negative) residual values, that is of importance in the test of preconditioning.

Figure 10 .
Figure 10.Matrix of maps showing potential predictor variables used in logistic regression analysis of hillslope failure probability.Each map shows distributed values of each predictor variable across the 5629 km 2 combined area of landslide mapping for both events (as shown in Fig.1).Variable descriptors and units are summarised in Table3.

Figure 12 .
Figure 12.Relative contributions of predictor variables to the fit of the 1929 and 1968 hillslope failure probability models.Sequence of model input predictors and resulting pseudo-R 2 goodness of fit values, produced by sequentially removing the least contributing predictor variable.(a) Results from Eq. (5) for the 1929 earthquake, (b) results from Eq. (5) for the 1968 earthquake, and (c) results from Eq. (6) for the 1968 earthquake.

Figure 13 .
Figure 13.Distributions of P LS residuals for the 1968 earthquake, hindcast using Eq.(6).Panels (a)-(c) show the residuals for this model (observed P LS minus predicted P LS ) plotted against each of the model predictors.Panel (d) shows the model residuals plotted against the predicted P LS .Panels (e) and (f) show the model residuals plotted against predicted P LS for the 1929 earthquake and distance from the 1929 coseismic fault, respectively.All residuals are calculated by aggregating probabilities across 20 equal quantile bins of the x variable.Positive residuals indicate that the model under-predicts P LS and negative residuals indicate that the model over-predicts P LS .Note that the amplitude of the plotted residuals varies due to binning and aggregating probabilities with the different predictor variables.For each plot, the coefficient (r) and significance (p) of correlation in the residuals are given.These were derived by adding each x variable into a logistic regression analysis of P LS predicted using Eq.(6) and observed landsliding.In this respect we test the coupled significance of SL and G as they feature in the model.

Figure C1 .
Figure C1.Matrix of variance inflation factors (VIF) for pool of all potential predictor variables.Combinations for which VIF values exceed the threshold of 10 (Kutner et al., 2004) are highlighted.
Name Date Epicentre location Magnitude Focal depth Rupture length Strike Dip Dip direction Buller 17 June 1929 41.70 • S, M s = 7.

Table 2 .
Summary of proxy variables suggested to influence spatial distributions of earthquake-induced landslides, based on empirical studies.
Seismic wave attributes (e.g.PGA, PGV, PGD, Local metric of shaking intensity

Table 3 .
Potential predictor variables, ID codes, descriptions and units.
(Kutner et al., 2004)ar coefficient of determination of the relationship between any two predictor variables.VIF values greater than 10 indicate a high level of multicollinearity and are avoided in our model(Kutner et al., 2004).The matrix of VIF values is given in Appendix C, and it indicates no high multicollinearity among the variables that feature in our final models.

Table 4 .
Logistic regression output coefficients and fit statistics.

Table 5 .
Logistic regression output coefficients and fit statistics, for models including the influence of the 1929 earthquake on the 1968 landslide distribution.