Odin Marc
email: odin.marc@get.omp.eu

Christoph Sens-Schönfelder

Luc Illien

Patrick Meunier

Manuel Hobiger

Kaoru Sawazaki

Claire Rault

Niels Hovius

Towards quantifying landscape mechanical variations after earthquakes with seismic interferometry

In mountainous terrain, large earthquakes often cause widespread coseismic landsliding as well as hydrologic disturbances. A subsequent transient phase with high landslide rates has also been reported for several earthquakes. Separately, subsurface seismic velocities are frequently observed to drop coseismically and subsequently recover. Consistently with various laboratory work we hypothesize that the seismic velocity changes may track coseismic damage and progressive recovery of landscape substrate, which modulate landslide hazard and hydrological processes, on timescales of months to years. To test this, we analyze the near-surface seismic velocity variations, obtained with single-station high-frequency (0.5 -4 Hz) passive image interferometry, in the epicentral zones of four shallow earthquakes, for which constraints on landslide susceptibility through time exist. In the case of the 1999 Chi-Chi earthquake, detailed landslide mapping allows us to accurately constrain an exponential recovery of landslide susceptibility with a relaxation timescale of about one year, closely matching the recovery of seismic velocities. The 2004 Niigata, 2008 Iwate, and 2015 Gorkha earthquakes have less resolved constraints on landsliding, but assuming an exponential recovery, we also find matching relaxation timescales, from ~0.1 to ~0.6 yr, for the landslide and seismic recoveries.

These observations support our hypothesis and suggest that systematic monitoring of seismic velocities after large earthquakes may help constrain and manage the evolution of landslide hazard in epicentral areas. To achieve this goal, we end by discussing several ways to tighten the link between seismic velocity and landscape mechanical properties, such as accounting for station-dependent velocity sensitivity to tidal strain, or a comparison with temporal evolution of permeability and/or debris-flow threshold as additional marker of landscape properties.

Introduction

Earthquakes greater than M w 5 can cause widespread landsliding [START_REF] Harp | Anomalous Concentrations of Seismically Triggered Rock Falls in Pacoima Canyon: Are They Caused by Highly Susceptible Slopes or Local Amplification of Seismic Shaking[END_REF][START_REF] Keefer | Investigating Landslides Caused by Earthquakes -A Historical Review[END_REF][START_REF] Khazai | Evaluation of factors controlling earthquake-induced landslides caused by Chi-Chi earthquake and comparison with the Northridge and Loma Prieta events[END_REF][START_REF] Meunier | Regional patterns of earthquake-triggered landslides and their relation to ground motion[END_REF][START_REF] Marc | A seismologically consistent expression for the total area and volume of earthquake-triggered landsliding[END_REF]Tanyas et al., 2017] as well as pervasive hydrological disturbance [START_REF] Wang | Field relations among coseismic ground motion, water level change and liquefaction for the 1999 Chi-Chi (Mw = 7.5) earthquake, Taiwan[END_REF][START_REF] Wang | Coseismic release of water from mountains: Evidence from the 1999 (Mw = 7.5) Chi-Chi, Taiwan, earthquake[END_REF][START_REF] Wang | Large earthquakes create vertical permeability by breaching aquitards[END_REF][START_REF] Wang | Hydrologic Responses to Earthquakes and a General Metric[END_REF]. These effects are mostly considered as an instantaneous response to static and dynamic deformation at the time of the earthquake. However, in the last decade, studies of groundwater well levels have shown that earthquakes can cause a long-lasting increase in permeability of the shallow subsurface (<100m) that recovers to pre-seismic levels over months to years [START_REF] Elkhoury | Seismic waves increase permeability[END_REF][START_REF] Xue | Continuous Permeability Measurements Record Healing Inside the Wenchuan Earthquake Fault Zone[END_REF][START_REF] Shi | Mechanism of co-seismic water level change following four great earthquakesinsights from co-seismic responses throughout the Chinese mainland[END_REF][START_REF] Shi | Aquifers switched from confined to semiconfined by earthquakes[END_REF]Wang et al., 2016). Similarly, a transient increase of landslide rates (normalized for the meteorological forcing) and recovery over similar timescales has been found for several large, shallow earthquakes [START_REF] Marc | Transient changes of landslide rates after earthquakes[END_REF][START_REF] Marc | Long-term erosion of the Nepal Himalayas by bedrock landsliding: the role of monsoons, earthquakes and giant landslides[END_REF]. These observations indicate that hydromechanical properties of landscapes can change substantially during and after earthquakes, with important implications for hydrology, geomorphology and natural hazards. However, the processes driving these changes, their spatial distribution and temporal evolution remain poorly constrained.

Passive seismic monitoring may offer a tool for characterization and monitoring of hydromechanical perturbations of landscape substrates, exploiting the dependency of seismic waves on the mechanical properties of the medium through which they travel. In fact, seismic interferometry based on ambient noise can retrieve variations of bulk seismic velocity (V) with high precision, dV/V ~0.1%, and at sub-daily temporal resolution. In the last decade, this approach has revealed and quantified a wide range of phenomena affecting the shallow crust down to several kilometers of depth, ranging from thermal stresses to water content (Snieder and[START_REF] Obermann | Depth sensitivity of seismic coda waves to velocity perturbations in an elastic heterogeneous medium[END_REF]Brenguier, 2019). Notably, many studies have reported seismic velocity drops during earthquakes, consistently followed by a progressive, non-linear recovery over months to years [START_REF] Brenguier | Postseismic Relaxation Along the San Andreas Fault at Parkfield from Continuous Seismological Observations[END_REF][START_REF] Wu | Temporal Changes in Site Response Associated with the Strong Ground Motion of the 2004 Mw 6.6 Mid-Niigata Earthquake Sequences in Japan[END_REF][START_REF] Sens-Schönfelder | Passive image interferometry for monitoring crustal changes with ambient seismic noise[END_REF][START_REF] Takagi | Coseismic velocity change in and around the focal region of the 2008 Iwate-Miyagi Nairiku earthquake[END_REF][START_REF] Wu | Long-term change of site response after the M W 9.0 Tohoku earthquake in Japan[END_REF][START_REF] Sawazaki | Time-lapse changes of P-and S-wave velocities and shear wave splitting in the first year after the 2011 Tohoku earthquake, Japan: shallow subsurface[END_REF][START_REF] Hobiger | Single-station cross-correlation analysis of ambient seismic noise: application to stations in the surroundings of the 2008 Iwate-Miyagi Nairiku earthquake[END_REF][START_REF] Hobiger | Coseismic and postseismic velocity changes detected by Passive Image Interferometry: comparison of one great and five strong earthquakes in Japan[END_REF][START_REF] Richter | Comprehensive observation and modeling of earthquake and temperature-related seismic velocity changes in northern Chile with passive image interferometry[END_REF][START_REF] Gassenmeier | Field observations of seismic velocity changes caused by shaking-induced damage and healing due to mesoscopic nonlinearity[END_REF]. Often, this cannot be attributed to changes of the static stress state or ground water level because these would generally cause the seismic velocity to decrease in some part of the epicentral area and increase in others, the latter not being observed [ e.g., Wegler et al. 2009, Sawazaki and[START_REF] Sawazaki | Time-lapse changes of P-and S-wave velocities and shear wave splitting in the first year after the 2011 Tohoku earthquake, Japan: shallow subsurface[END_REF][START_REF] Hobiger | Single-station cross-correlation analysis of ambient seismic noise: application to stations in the surroundings of the 2008 Iwate-Miyagi Nairiku earthquake[END_REF][START_REF] Hobiger | Coseismic and postseismic velocity changes detected by Passive Image Interferometry: comparison of one great and five strong earthquakes in Japan[END_REF]. Instead, wide-spread shallow substrate damage and subsequent healing have been expressed as a possible cause [e.g., [START_REF] Rubinstein | Depth constraints on nonlinear strong ground motion from the 2004 Parkfield earthquake[END_REF][START_REF] Sleep | Depth of Rock Damage from Strong Seismic Ground Motions near the 2004 Parkfield Mainshock[END_REF][START_REF] Takagi | Coseismic velocity change in and around the focal region of the 2008 Iwate-Miyagi Nairiku earthquake[END_REF][START_REF] Hobiger | Single-station cross-correlation analysis of ambient seismic noise: application to stations in the surroundings of the 2008 Iwate-Miyagi Nairiku earthquake[END_REF][START_REF] Hobiger | Coseismic and postseismic velocity changes detected by Passive Image Interferometry: comparison of one great and five strong earthquakes in Japan[END_REF][START_REF] Sawazaki | Time-lapse changes of P-and S-wave velocities and shear wave splitting in the first year after the 2011 Tohoku earthquake, Japan: shallow subsurface[END_REF][START_REF] Richter | Comprehensive observation and modeling of earthquake and temperature-related seismic velocity changes in northern Chile with passive image interferometry[END_REF][START_REF] Gassenmeier | Field observations of seismic velocity changes caused by shaking-induced damage and healing due to mesoscopic nonlinearity[END_REF] Several aspects prevent a straightforward attribution of any observed seismic velocity change to landscape mechanical properties. First, given that seismic interferometry retrieves the mean seismic velocity in a volume of crust, it is difficult to determine if changes occur mainly in the near-surface (<100-200m) or at depth, or uniformly throughout the crust. This leaves room for radically different interpretations of the same data (e.g., [START_REF] Brenguier | Postseismic Relaxation Along the San Andreas Fault at Parkfield from Continuous Seismological Observations[END_REF][START_REF] Sleep | Depth of Rock Damage from Strong Seismic Ground Motions near the 2004 Parkfield Mainshock[END_REF]. For example, velocity drop and recovery have alternatively been attributed to crustal stress redistribution [START_REF] Brenguier | Postseismic Relaxation Along the San Andreas Fault at Parkfield from Continuous Seismological Observations[END_REF][START_REF] Wegler | Sudden drop of seismic velocity after the 2004 Mw 6.6 mid-Niigata earthquake, Japan, observed with Passive Image Interferometry[END_REF][START_REF] Obermann | Seismic noise correlations to image structural and mechanical changes associated with the Mw 7.9 2008 Wenchuan earthquake[END_REF], 2018[START_REF] Wang | Evidence of Changes of Seismic Properties in the Entire Crust Beneath Japan After the Mw 9.0, 2011 Tohoku-oki Earthquake[END_REF] or weakening of fluid-pressurized parts of the crust [START_REF] Brenguier | Mapping pressurized volcanic fluids from induced crustal seismic velocity drops[END_REF]. Second, the physical mechanism underlying what is referred to as damage or weakening of rocks is expected to vary with material properties (e.g., granular material, or rock) and with the physical conditions related to depth (pressure, temperature). Laboratory experiments have illustrated a large range of non-linear response and recovery related to moderate strain or shaking of rocks and granular materials [START_REF] Johnson | Resonance and elastic nonlinear phenomena in rock[END_REF][START_REF] Ten Cate | Slow dynamics in the nonlinear elastic response of Berea sandstone[END_REF][START_REF] Tencate | Universal Slow Dynamics in Granular Solids[END_REF][START_REF] Tremblay | Probing slow dynamics of consolidated granular multicomposite materials by diffuse acoustic wave spectroscopy[END_REF][START_REF] Jia | Elastic weakening of a dense granular pack by acoustic fluidization: Slipping, compaction, and aging[END_REF][START_REF] Rivière | A set of measures for the systematic classification of the nonlinear elastic behavior of disparate rocks[END_REF]. Notably, equivalent results were obtained during field experiments of ground shaking [START_REF] Lawrence | Situ Measurements of Nonlinear and Nonequilibrium Dynamics in Shallow, Unconsolidated Sediments[END_REF][START_REF] Renaud | In situ characterization of shallow elastic nonlinear parameters with Dynamic Acoustoelastic Testing[END_REF]. These observations have been interpreted with various physical models incorporating non-linear elasticity of grain-to-grain or fracture contacts, [START_REF] Guyer | Nonlinear Mesoscopic Elasticity: Evidence for a New Class of Materials[END_REF][START_REF] Broda | Modelling of nonlinear crackwave interactions for damage detection based on ultrasound-A review[END_REF][START_REF] Snieder | The time dependence of rock healing as a universal relaxation process, a tutorial[END_REF], Sens-Schönfelder et al., 2019). Similarly, material damage and recovery observed at larger strains and confining pressures have been interpreted in terms of fracture mechanics [START_REF] Brantut | Time-dependent recovery of microcrack damage and seismic wave speeds in deformed limestone[END_REF]. Third, environmental changes, including temperature, precipitation, groundwater level fluctuation, and freeze-and-thaw can cause substantial, often quasi-cyclic, variations of seismic velocity [START_REF] Sens-Schönfelder | Passive image interferometry and seasonal variations of seismic velocities at Merapi Volcano, Indonesia[END_REF], Hillers et al., 2014[START_REF] Richter | Comprehensive observation and modeling of earthquake and temperature-related seismic velocity changes in northern Chile with passive image interferometry[END_REF][START_REF] Wang | Seasonal Crustal Seismic Velocity Changes Throughout Japan[END_REF][START_REF] Clements | Tracking Groundwater Levels Using the Ambient Seismic Field[END_REF][START_REF] Fores | Monitoring Saturation Changes with Ambient Seismic Noise and Gravimetry in a Karst Environment[END_REF] that can overprint any changes associated with earthquakes [START_REF] Hobiger | Coseismic and postseismic velocity changes detected by Passive Image Interferometry: comparison of one great and five strong earthquakes in Japan[END_REF][START_REF] Wang | Seasonal Crustal Seismic Velocity Changes Throughout Japan[END_REF]. At the scale of a single landslide, recent observations suggested important coupling between shaking and rainfall events modulated the recovery of the materials [START_REF] Bontemps | Rain and small earthquakes maintain a slowmoving landslide in a persistent critical state[END_REF].

For these reasons, the ability to monitor mechanical landscape properties relevant to surface processes with seismological means has remained unclear, and a formal comparison of observed velocity changes with independent estimates of landscape mechanical property changes has not yet been documented. Here we compare new and published seismological and geomorphic data from four earthquakes, to test the hypothesis that velocity changes measured by seismic interferometry during and after large, shallow earthquakes, can be attributed to changes in the hydromechanical properties of landscape substrates. This hypothesis is plausible given the laboratory and field experiments reviewed above. It can be supported or rejected based on the correspondence between the seismic and geomorphic recovery timescales and behavior, including the magnitude of the post-earthquake variations of near-surface seismic velocities and of independent proxies for landscape properties such as landslide susceptibility.

We show that, within their limitations, the existing data yield consistent recovery timescales, and highlight a well-documented case where this recovery followed, an exponential trajectory with a measurable decay constant. For two cases, available depth constraints are consistent with maximal velocity changes near the surface (<~100m). These findings constitute preliminary confirmation of our hypothesis. We end by discussing future (re-)analyses that would allow us to better test our hypothesis.

STUDY AREAS, DATASET AND GENERAL METHODS

Comparing mechanical properties of the subsurface that affect landslide occurrence with elastic properties sensed by seismic waves, imposes severe constraints on the dataset selection. It requires a seismically instrumented epicentral area of a shallow earthquake in mountainous terrain with sufficient landsliding and satellite visibility. Time-dependent changes in landslide susceptibility have been observed and quantified in the epicentral areas of several large, shallow earthquakes including the 1993 Mw 6.9 Finisterre earthquake in Papua New Guinea, 1999 Mw 7.6 Chi-Chi earthquake in Taiwan, the 2004 Mw 6.6 Mid-Niigata and the 2008 Mw 6.9 Iwate-Miyagi earthquake in Japan [START_REF] Marc | Transient changes of landslide rates after earthquakes[END_REF] and the 2015 Mw 7.9 Gorkha earthquake in Nepal [START_REF] Marc | Long-term erosion of the Nepal Himalayas by bedrock landsliding: the role of monsoons, earthquakes and giant landslides[END_REF]. In all these areas except Papua New Guinea, available local instrument records allow the assessment of seismic velocity variations by passive image interferometry (Sens-Schönfelder and Wegler, 2011, Fig 1-4), based on the correlation of the ambient seismic wavefield. However, of all these cases, only the Chi-Chi earthquake has a highly resolved temporal evolution of landslide susceptibility, that can be compared with the evolution of shallow seismic velocities. Therefore, after introducing the general methodology, we first present detailed results of this exceptional case, and use it as a benchmark to analyze the observations of the three other cases (Niigata, Iwate and Gokha earthquakes) in a second time.

In order to resolve and compare the transient responses of geophysical parameters after an earthquake, we focus on continuous data, with high temporal resolution spanning several years to a decade. Furthermore, we aimed at describing the representative behavior of the epicentral areas, by averaging the responses of several seismic stations and the regional occurrence of landslides. Note that the spatial extent is different between the data types used in this study, and that any comparison of time-dependent behavior assumes that the processes affecting the materials near the Earth's surface are similar for the areas that are sampled by high-frequency (0.5-4Hz) seismic waves recorded at a seismometer and the hillslopes affected by distributed landsliding.

Landslide susceptibility changes

In the steep mountain topography of our study areas, widespread landsliding can occur during strong rainfall events (typhoons or monsoons) or earthquakes. At other times, the landslide rate is much lower. Here we use published landslide inventories [START_REF] Marc | Transient changes of landslide rates after earthquakes[END_REF], 2019, Fig 2, 5, S1), containing all new landslides that occurred within time intervals constrained by the availability of cloud-free imagery, spanning a few weeks to a few months. During these intervals, conditions ranged from sparse rainfall to extreme typhoon precipitation. This variability of rainfall forcing on landslide triggering must be taken into account to obtain a metric characterizing the mechanical state of the landscape subsurface. We define the landslide susceptibility as the amount of landsliding (volume or area) in a time interval normalized by the expected amount of landsliding for the magnitude of the rainfall during the interval. This ratio requires calibration of a local relationship between landsliding and rainfall. In Japan and Taiwan, consistently with the assumption that rainfall is the main landslide trigger, an exponential function of the total rainfall during the mapping interval adequately explains the landsliding observed during all the intervals preceding the mainshock, and the ones occuring several years after (Marc et al., 2015, Fig 2, S1). In contrast, periods immediately following an earthquake systematically had much more landsliding than other intervals with similar rainfall forcing [START_REF] Marc | Transient changes of landslide rates after earthquakes[END_REF]. We note that shallow, large (Mw >5) aftershocks can be an additional trigger of landsliding, but in the cases considered here, all large aftershocks occurred within a period shorter than the episode of elevated post-seismic landslide susceptibility [START_REF] Marc | Transient changes of landslide rates after earthquakes[END_REF].

In Nepal, landslides are often due to sustained monsoon rain, rather than shorter precipitation events. [START_REF] Marc | Long-term erosion of the Nepal Himalayas by bedrock landsliding: the role of monsoons, earthquakes and giant landslides[END_REF] mapped landslides over an eight-year period in four valleys, three substantially affected by the Gorkha earthquake rupture, with many coseismic landslides, and one, the Kali Gandaki valley, unaffected by this earthquake as it is >100 km west of the rupture plane and opposite to directivity effect ( Fig 1, Roback et al, 2017). They pointed out three deep-seated landslides distributed over two valleys, which lie outside the normal sizedistribution statistics and were likely due to long-term gravity-driven slope deformation, rather than any increase in pore-pressure due to rainfall (see Fig 4 in [START_REF] Marc | Long-term erosion of the Nepal Himalayas by bedrock landsliding: the role of monsoons, earthquakes and giant landslides[END_REF]. They also excluded many large slope failures caused by a glacial lake outburst flood in one of the valleys in 2016 (Cook et al., 2018), as well as the ones due to the largest aftershock (Marc et al., 2019, Fig 1). After these corrections, estimated landslide rates were found to be relatively steady, and similar across the four different valleys, but a relationship with satellite-based local monsoon rainfall did not emerge. Therefore, we used as a base line the landsliding in the Kali Gandaki valley, which was not affected by strong ground motion during the Gorkha earthquake, but is otherwise comparable to the three other valleys directly above the part of the fault that slipped in 2015 [START_REF] Marc | Long-term erosion of the Nepal Himalayas by bedrock landsliding: the role of monsoons, earthquakes and giant landslides[END_REF].

Seismic velocity changes

Relative seismic velocity variations, dV/V, were observed with passive image interferometry [START_REF] Sens-Schönfelder | Passive image interferometry for monitoring crustal changes with ambient seismic noise[END_REF], based on the correlation of the ambient seismic wavefield (Fig. 234). This technique allows for the continuous, high-resolution (~daily) retrieval of the average dV/V in a volume of rocks of kilometric scale surrounding the stations.

For the documented earthquakes in Taiwan and Japan, we performed analyses in the 0.5-1, 1-2 and 2-4Hz frequency bands, using single-station cross-correlation [with data processing identical to [START_REF] Hobiger | Single-station cross-correlation analysis of ambient seismic noise: application to stations in the surroundings of the 2008 Iwate-Miyagi Nairiku earthquake[END_REF], because of the superior sensitivity of variations in the high frequencies at shallow depths. For the Gorkha earthquake, Nepal, we had access to data from a one component station (Z) allowing auto-correlation only. The processing was identical to the other stations except that we used a coda time range of 2-14 times the longest period of each frequency range, instead of 7.5 -17.5 times as in the Hobiger method. This was done because the energy levels of the later records of this station were very low. This difference in processing detail should not have a major bearing on the reported results.

The coseismic and post-seismic dV/V in the Iwate and Niigata epicentral areas have been described earlier [START_REF] Takagi | Coseismic velocity change in and around the focal region of the 2008 Iwate-Miyagi Nairiku earthquake[END_REF][START_REF] Hobiger | Single-station cross-correlation analysis of ambient seismic noise: application to stations in the surroundings of the 2008 Iwate-Miyagi Nairiku earthquake[END_REF][START_REF] Hobiger | Coseismic and postseismic velocity changes detected by Passive Image Interferometry: comparison of one great and five strong earthquakes in Japan[END_REF]. Here, we present our analysis of the velocity changes at stations ICEH / ICWH and YNTH / NGOH in the focal zone of the Iwate-Miyagi and Niigata earthquakes, respectively. They allow robust retrieval of strong velocity variations across several frequency bands, showing an increase of the magnitude of the coseismic dV/V drop with increasing frequency (Fig 3), consistent with large changes at shallow depth(~5%, <~200m) [START_REF] Takagi | Coseismic velocity change in and around the focal region of the 2008 Iwate-Miyagi Nairiku earthquake[END_REF][START_REF] Hobiger | Single-station cross-correlation analysis of ambient seismic noise: application to stations in the surroundings of the 2008 Iwate-Miyagi Nairiku earthquake[END_REF]. We do not consider station NRKH, which also recorded important velocity changes, but with a similar magnitude across frequencies. This has been shown to correspond to diffuse changes over a greater depth range (~0.1%, up to ~3-4 km) [START_REF] Takagi | Coseismic velocity change in and around the focal region of the 2008 Iwate-Miyagi Nairiku earthquake[END_REF][START_REF] Hobiger | Single-station cross-correlation analysis of ambient seismic noise: application to stations in the surroundings of the 2008 Iwate-Miyagi Nairiku earthquake[END_REF], possibly near the fault plane at depth. For the Chi-Chi earthquake in Taiwan we have used records from seismic stations SSLB and TDCB, the only stations in the Broadband Array in Taiwan for Seismology (BATS) array that were located less than 50km away from the rupture plane recording before the main shock.

The TDCB signal did not exhibit significant velocity changes after the earthquake (Fig S2). For SSLB, the stability of noise correlation functions permitted retrieval of velocity changes within the 2-4Hz range and for the horizontal components (E-E, E-N, N-N) only. In Nepal, we could only access ~2 years of records from 01/01/2014 to 09/12/2015 from the GUMBA station, managed by the Nepal Seimological Center. This station is located on a ridge, within the Bhote Koshi valley, directly above the rupture plane of the Gorkha earthquake. We could not retrieve velocity changes at 0.5-1Hz but retrieved stable correlations and substantial coseismic velocity changes in the 1-2 and2-4Hz frequency bands (Fig 4).

All selected stations (SSLB, ICWH, ICEH, YNTH, NGOH and GUMBA) are within or very close to an area where landslides triggered by the respective earthquakes were mapped ( Fig 1).

Determination of relaxation times and confidence intervals

Comparison of relaxation times is an important aspect of our study. Here we detail the approaches used to constrain the time scale of the transient response of the different observables, at the scale of the epicentral area.

For seismic velocities, we used models containing an exponential post-earthquake relaxation, well-suited to describe the long-term mechanical recovery [START_REF] Snieder | The time dependence of rock healing as a universal relaxation process, a tutorial[END_REF], with a seasonal cycle (with annual pulsation w=2pi yr -1 ), and a residual change:

where O is the observable, A the initial value, S and  the amplitude and phase of the seasonal signal, C the permanent coseismic change, P the recovering change, τ the relaxation time, t E the time of the earthquake and H the Heaviside function. We fit the landslide susceptibility with a similar function, imposing A=1, and S=0 because the susceptibility measure already includes a normalization for (seasonal) rainfall variations. The uncertainty of any constraints on the background landslide susceptibility does not allow an accurate assessment of permanent (longterm) changes in landslide susceptibility. However, such changes are likely smaller than the typical uncertainties arising from the rainfall normalization (about a factor of two), and we assume C=0, meaning the landslide susceptibility only depends on P and τ and returns to one after a sufficient amount of time. Below, we use the subscripts S and L on P and τ to refer to the seismic velocity or landslide susceptibility, respectively.

Seismic relaxation time and their confidence intervals were obtained by using a nonlinear least-square fitting of models with given relaxation time varied logarithmically from 0.01 to 10 years. For each model variance of the residuals and variance ratio relative to the best model was computed. The variance ratio at which models are statistically indistinguishable from the best model at 95% confidence level (which increases with the amount of data used to constrained the model) were assessed with an F-test, allowing us to define confidence interval around the best relaxation time. To obtain one relaxation time estimate and interval at each seismic station with results over several frequency bands, we merged residuals from all frequencies and performed an F-Test on those (Fig. 3,4). Then, in Japan, we consider the mean of the two best estimate and intervals as the most representative values for the epicentral area.

TIME-DEPENDENT EVOLUTION OF GEOMORPHIC AND SEISMIC CHANGES

The detailed case of the Chi-Chi earthquake

The Chi-Chi earthquake increased the regional landslide susceptibility by a factor larger than 20. We modeled the landslide susceptibility with Eq 1 (with A=1, S=C=0) and fitted the mean model prediction in a time series consisting of 32 mapping intervals to the observed mean landslide susceptibility. To estimate uncertainties on P L and  S , we repeated this fitting procedure 1000 times, while randomly selecting P L between 20 and 40 and multiplying each susceptibility estimate by a random coefficient between 0.5 and 2, to account for the typical uncertainties.

With this procedure we obtained a post-seismic change magnitude P L between 28 and 32 and a median relaxation time constant τ L of 0.9 yr, (0.8-1, the 5th-95th percentiles of the modeled relaxation) (Fig. 2). During this recovery, thousands of new landslides occurred, mostly on pristine slopes, and mostly due to rainfall rather than further seismic activity [START_REF] Marc | Transient changes of landslide rates after earthquakes[END_REF].

Meanwhile, noise correlation at station SSLB, at the northern edge of the area where landslides were mapped, revealed a dV/V of about -1% in the 2-4Hz band, followed by a non-linear recovery (Fig. 2). Lacking a local subsurface seismic velocity model, we assume that the station has a depth sensitivity in the 2-4Hz band similar to that of the Japanese stations (<0.2-2km, see 3.2). This depth range is a first order estimate only, but it comprises the uppermost 10 m of the solid Earth, in which most post-seismic landslides root, as estimated from scaling relationships between landslide area and depth (Larsen et al., 2010) [cf., Fig 3 of Marc et al., 2015]. This depth range excludes a predominance of effects of near-fault changes. A strong seasonal signal overprints the post-seismic dV/V trend, likely related to seasonal, cyclic variations of the groundwater level [START_REF] Sens-Schönfelder | Passive image interferometry and seasonal variations of seismic velocities at Merapi Volcano, Indonesia[END_REF]Hillers et al., 2014]. We varied τ from 0.1 to 10 years while optimizing the other fit parameters to obtain a best fit time constant for the recovery of the dV/V of τ S =1. 1 [0.45-3.1] yr (F-test 95% Confidence interval, Fig 2,7).

The Chi-Chi earthquake case is sufficiently well documented to constrain both seismic and landslide behavior with confidence. Importantly the landslide time-series strongly suggests an exponential relaxation with a relaxation time matching the best estimates obtained for the dV/V recovery. We consider this case as a benchmark, and for the other cases with more limited data we test whether a similar exponential model applied for both time-series also yields similar timescales. Indeed, in the three other cases discussed here, the constraints on the recovery of landslide susceptibility are considerably less precise, primarily because the shorter relaxation periods did not permit to acquire landslide data in a larger number of time intervals.

However, for these cases, good seismic data exist, motivating us to consider them here.

Complimentary constraints from the Niigata, Iwate and Gorkha earthquakes

Constraining landscape recovery with few post-seismic intervals

For the Niigata, Iwate and Gorkha cases, only one landslide mapping intervals after the earthquake exhibited an increased landslide susceptibility. This does not allow a robust estimation of the recovery timescale. Instead, assuming an exponential recovery similar to the one of the Chi-Chi case, we explored the range of τ L that can fit the average landslide susceptibility during both the first and second post-seismic mapping intervals (Fig 6). P L is under-constrained, and large values would require very small relaxation times to match the data. However, it is reasonable to expect that for the Iwate-Miyagi, Niigata and Gorkha cases, the initial perturbation must have been of the same order of magnitude or smaller than that of the Chi-Chi earthquake, which had stronger ground shaking and P L of around 30 (Fig 2A). We therefore varied P L between 10 and 30 (Fig 6). Accounting for a factor of 2 uncertainty on the landslide susceptibility [START_REF] Marc | Transient changes of landslide rates after earthquakes[END_REF] we obtain a broader range of values of  L that bracket what we consider to be the most likely range (when considering the landslide susceptibility at face value). For these cases we refrain from giving a best estimate.

Two Japanese cases with constraints on landslides

At high frequencies (2-4Hz), dV/V dropped by 0.8-1% around stations in the epicentral areas of the Niigata and Iwate earthquakes in Japan, and the average relaxation times of the dV/V were τ S =0.32 yr [0.16-0.49, F-test 95% Confidence interval] for the former, and τ S =0.83 [0.52-1.6] yr for the latter (Fig 3 ,7). Notably, seismic observations from distributed stations for these two cases show that dV/V relaxation times varied somewhat (Fig 3) within the epicentral area. [START_REF] Hobiger | Single-station cross-correlation analysis of ambient seismic noise: application to stations in the surroundings of the 2008 Iwate-Miyagi Nairiku earthquake[END_REF] computed dispersion curves based on perturbed seismic velocity models at different frequencies to retrieve the approximate volumes within which these velocity changes occurred. They found that, for stations ICEH and ICWH, velocity drops of 5% between 0 -100 / 200 m depth produced a good fit of the data, as did smaller drops of 0.3-1% extending down to 1-2 km depth. Using this approach for station YNTH, which is the only one with a local velocity model, our best estimate is a velocity change of 0.5% between the surface and 1.5-3 km depth, although a 2% drop in the first 100m has an almost equally low misfit (Fig S3).

For the Niigata and Iwate earthquakes, the models best matching the observed landslide susceptibility have τ L = 0.08-0.25 yr [0.05-0.3 yr, considering a factor two uncertainty on the susceptibilities] with P L unconstrained between 10 and 30, and τ L = 0.3-0.5 [0.15-0.7] yr with P L likely between 20 and 30, respectively. Thus, for both Japanese cases, τ S overlapped with the corresponding τ L, although it is unclear for now if the slightly faster geomorphic relaxation in these two cases is physical or an artifact of limited data.

The Gorkha earthquake case

We found that the mean total landsliding (in area and volume) in the three valleys affected by the earthquake was broadly proportional to the landsliding in the Kali Gandaki valley, except during the interval just after (but not containing) the mainshock, in 2015 (Fig 5A).

Hence, we define a mean landslide susceptibility across the epicentral area, averaging the landsliding in the three affected valleys and normalizing with respect to the mapped landsliding in the Kali Gandaki valley (Fig 5B). This normalization defines a strong, 18-fold landslide susceptibility increase during the 2015 summer monsoon season, shortly after the Gorkha earthquake, and a rapid decay to near background susceptibility in the next monsoon season.

Thus, the geomorphic recovery in the Nepali case is more strictly constrained than in Japan, assuming P L is not much larger than 30. Thus we obtained τ L = 0.25-0.35 [0.18-0.42] yr, with P L most likely near 30.

For this case, we observed a strong dV/V drop of 6% at 2-4 Hz, followed by a partial recovery of ~45% of this change within 100 days, without any evidence of substantial seasonal cycles (Fig 4). At 1-2Hz, the dV/V drop was smaller, about 1.5%, and in this frequency range substantial seasonal velocity variations were evident in both years. In this frequency range, around 70% of the coseismic change was recovered by the end of the record, about 210 days after the earthquake. The short span of the data does not allow us to constrain with high confidence whether the Gorkha earthquake caused a permanent offset of the seismic velocity, nor if it is meaningful to compare between the seasonal cycles of 2014 and 2015, even though they appear quite different at face value. We note that the 1-2 Hz dV/V declined steadily for almost two months after the earthquake, possibly due to effects of aftershocks. This interval included an Mw 7.3 event, 17 days after the mainshock, with a source located less than 30 km away from station GUMBA. Explaining or correcting for these effects is beyond the scope of this study, but our observations suggest that longer data coverage, and a more explicit modeling of the impact of rainfall and aftershocks on the dV/V could refine the constraints on the relaxation time. For now, applying our simple methods as for other cases, we obtained a relaxation time τ S =0.16 yr [0.08-0.5 yr, F-test 95% Confidence interval], mainly driven by the longer, more complete recovery of the 1-2Hz result (Fig 4C). In contrast to the three other cases, the seismic relaxation appears to be somewhat shorter than the landslide relaxation time, but the estimated ranges overlap (Fig 7).

DISCUSSION

Depth constraints on the subsurface velocity change

Landslide susceptibility of steep landscapes is substantially controlled by material properties within around 10m below the surface, where almost all landslides occur (Larsen et al., 2010). It is also set by substrate permeability changes, which are often reported within 100-200m from the surface (e.g., [START_REF] Elkhoury | Seismic waves increase permeability[END_REF][START_REF] Xue | Continuous Permeability Measurements Record Healing Inside the Wenchuan Earthquake Fault Zone[END_REF]. Of interest, therefore, are seismic velocity changes in the uppermost part of the subsurface. These changes should be distinguished from any effects on seismic velocity due to deep crustal processes, associated with deep fluids (e.g., [START_REF] Brenguier | Mapping pressurized volcanic fluids from induced crustal seismic velocity drops[END_REF][START_REF] Wang | Evidence of Changes of Seismic Properties in the Entire Crust Beneath Japan After the Mw 9.0, 2011 Tohoku-oki Earthquake[END_REF], near-fault damage or stress redistribution across large depth ranges [START_REF] Brenguier | Postseismic Relaxation Along the San Andreas Fault at Parkfield from Continuous Seismological Observations[END_REF], Qui et al., 2020). When velocity changes can be estimated for multiple frequency-bands, assuming or constraining the wavefield (surface waves or scattered waves) and a velocity model allows the retrieval of the most likely depth and amplitude of the change [START_REF] Obermann | Depth sensitivity of seismic coda waves to velocity perturbations in an elastic heterogeneous medium[END_REF][START_REF] Hobiger | Single-station cross-correlation analysis of ambient seismic noise: application to stations in the surroundings of the 2008 Iwate-Miyagi Nairiku earthquake[END_REF]Fig S3).

Although large trade-offs are involved, this approach allows the distinction between sites where dV/V is likely dominated by shallow changes (e.g., station ICEH) and other sites dominated by deep (>2km) changes (e.g., station NRKH located above the fault plane of the Iwate earthquake) [START_REF] Hobiger | Single-station cross-correlation analysis of ambient seismic noise: application to stations in the surroundings of the 2008 Iwate-Miyagi Nairiku earthquake[END_REF]. Perhaps the best opportunity for estimation of dV/V in a well-defined, shallow depth range is offered by KiK-net stations in Japan, which comprise a surface seismometer and a borehole seismometer, typically at 100-200m depth [START_REF] Takagi | Coseismic velocity change in and around the focal region of the 2008 Iwate-Miyagi Nairiku earthquake[END_REF][START_REF] Sawazaki | Time-lapse changes of P-and S-wave velocities and shear wave splitting in the first year after the 2011 Tohoku earthquake, Japan: shallow subsurface[END_REF][START_REF] Nakata | Estimating near-surface shear wave velocities in Japan by applying seismic interferometry to KiK-net data[END_REF]. At KiK-net stations, shallow dV/V values have been found to be much larger (>5%) than changes averaged over a greater depth range. Moreover, these sites allowed the comparison of analyses constrained by deconvolution within shallow boreholes and by noise correlation. Such studies, which include the stations we have re-analyzed for the Niigata and Iwate earthquakes, suggest that dV/V are mostly located within the first few hundred meters below the surface [START_REF] Sawazaki | Depth-dependence of seismic velocity change associated with the 2011 Tohoku earthquake, Japan, revealed from repeating earthquake analysis and finite-difference wave propagation simulation[END_REF][START_REF] Sawazaki | Estimation of seismic velocity changes at different depths associated with the 2014 Northern Nagano Prefecture earthquake, Japan (MW6.2) by joint interferometric analysis of NIED Hi-net and KiK-net records[END_REF]. For our two Japanese cases, we have found that the amplitude of the velocity drop increases systematically with the frequency of the seismic waves, while the post-seismic relaxation time does not (Fig 3). Taken together, these observations suggest that the dV/V we studied primarily result from changes in the very shallow layers, where shaking is most intense and confining pressure lowest. In view of this, velocity changes of 0.1-1%, obtained at low frequency (<1Hz) [e.g., [START_REF] Brenguier | Postseismic Relaxation Along the San Andreas Fault at Parkfield from Continuous Seismological Observations[END_REF][START_REF] Sens-Schönfelder | Passive image interferometry for monitoring crustal changes with ambient seismic noise[END_REF], may result from averaging over large depth ranges with small changes at large depths and changes of the order of 10% in the first tens of meters near the Earth's surface [START_REF] Sawazaki | Time-Lapse Changes of Seismic Velocity in the Shallow Ground Caused by Strong Ground Motion Shock of the 2000 Western-Tottori Earthquake, Japan, as Revealed from Coda Deconvolution Analysis[END_REF][START_REF] Takagi | Coseismic velocity change in and around the focal region of the 2008 Iwate-Miyagi Nairiku earthquake[END_REF][START_REF] Hobiger | Single-station cross-correlation analysis of ambient seismic noise: application to stations in the surroundings of the 2008 Iwate-Miyagi Nairiku earthquake[END_REF]. Such large changes in the shallow subsurface could be directly associated with measurable changes in the dynamics of hydrological or geomorphological processes. Therefore, we proceed to consider the physical processes that could explain landslide susceptibility changes in mountain landscapes induced by seismic strong ground motion.

Mechanical interpretation of landslide susceptibility increase and recovery

Seismically-induced increases in landslide susceptibility may be due to processes reducing the resisting forces, i.e., reduction of the material strength, or processes increasing the forces driving instability, i.e., pore-pressure increases or steepening of hillslopes. Loss of cohesion and/or friction, through seismically-induced ground cracking, and/or regolith decompaction, are often observed after earthquakes (Owen et al., 2008, Collins andJibson, 2015). Strength loss in surface materials may also result from vegetation death and root decay, as observed after tree-cutting (Sidle, 1991) or volcanic eruptions [START_REF] Fan | Earthquake-Induced Chains of Geologic Hazards: Patterns, Mechanisms, and Impacts[END_REF], but widespread decay of vegetation is not typically observed after earthquakes (Vittoz et al., 2001) and has been ruled out as a major factor. Furthermore, coseismic deformation may affect the subsurface hydrology but it does not necessarily favor higher pore pressures within the first few meters, where landslide failures usually locate. For example, it has been hypothesized that the probability of high pore-pressure in the shallow regolith is reduced by widespread cracking due to strong ground motion, and a resulting increase of the vertical permeability of near surface rocks, enhancing water drainage toward the groundwater table at depth [START_REF] Wang | Coseismic release of water from mountains: Evidence from the 1999 (Mw = 7.5) Chi-Chi, Taiwan, earthquake[END_REF][START_REF] Wang | Large earthquakes create vertical permeability by breaching aquitards[END_REF]. Finally, coseismic steepening of hillslopes could explain higher landslide susceptibility, but coseismic surface deformation does not substantially change the topography except in the immediate vicinity of fault surface ruptures and coseismic landslides. Most post-seismic landslides (70-90%) are not colocated with these features [START_REF] Marc | Transient changes of landslide rates after earthquakes[END_REF][START_REF] Marc | Long-term erosion of the Nepal Himalayas by bedrock landsliding: the role of monsoons, earthquakes and giant landslides[END_REF]. On balance, the coseismic increase of landslide susceptibility is most likely due to a spatially distributed reduction of the strength of near surface materials (due to various forms of rock and regolith damage). The magnitude of this strength reduction likely scales with P L . In addition, coseismic damage may also cause hydrological changes affecting landslide dynamics. Laboratory experiments confirm that such changes in geomechanical properties result in a substantial seismic velocity drop (e.g., [START_REF] Shokouhi | Dynamic stressing of naturally fractured rocks: on the relation between transient changes in permeability and elastic wave velocity[END_REF]. In the landscape, this effect is likely concentrated near the surface, as is evident from observations in Japan (see Depth constraints on the subsurface velocity change). An intuitive explanation of the observed post-seismic decrease (recovery) of the landslide susceptibility is the progressive removal of seismically damaged layers (Fig S4).

However, the influence of healing processes is implied by the recovery of the seismic velocity also in locations without notable erosion. This is the case at stations ICEH, NGOH and YNTH, that are located in valley floors where erosion is unlikely (Fig 1). It is likely that the healing of the shallow substrate that occurs at these stations also affects landslide-prone terrain, and it may be the dominant process driving the reduction and recovery of landslide susceptibility after earthquakes. Furthermore, we expect that the removal of a damaged layer from the surface of topography perturbed by seismic ground motion causes a reduction of the correlation (or coherence) of the seismic wavefield rather than a return to pre-seismic background velocities, given that it would change the scattering structure of the station surroundings (Planes et al., 2014[START_REF] Obermann | Seismic noise correlations to image structural and mechanical changes associated with the Mw 7.9 2008 Wenchuan earthquake[END_REF]. Therefore, a comparative analysis of the temporal evolution of the seismic wave velocity and the coherence of the waveform may help to further clarify the mechanisms explaining the recovery of landslide susceptibility. This analysis is beyond the scope of this paper.

Interpretation of quasi-permanent velocity changes

After the Chi-Chi earthquake, landslide susceptibility and seismic velocities appear to have co-evolved. In our other cases, these two variables also had similar recovery timescales, but they are less well constrained. However, we note that although the landslide rate appears to have recovered to background values in each case, the magnitude of the seismic velocity recovery varied between cases, from 30% to 100% of the coseismic drop, without a consistent trend with depth as inferred from frequency (Fig 2,3,6). The remaining discrepancy with the pre-earthquake seismic velocity persisted on the time scale of our observations, and we consider it quasi-permanent. Uncertainties in the normalization for the variable strength of rainfall as a trigger of landsliding may obscure a permanent change of landslide susceptibility, but they are typically between less than 10% (Chi-Chi, Gorkha) and 30% (Niigata) of the coseismic change. This is significantly smaller than the permanent change in seismic velocity observed in some of our cases (Fig 3,6).

A quasi-permanent velocity change may be due to plastic behavior when strain exceeds a threshold [START_REF] Brunet | Transitional nonlinear elastic behaviour in dense granular media[END_REF], and its variable magnitude may relate to the lithological and rheological diversity of the geomaterials in the epicentral areas and at different depths [START_REF] Takagi | Coseismic velocity change in and around the focal region of the 2008 Iwate-Miyagi Nairiku earthquake[END_REF][START_REF] Sawazaki | Time-lapse changes of P-and S-wave velocities and shear wave splitting in the first year after the 2011 Tohoku earthquake, Japan: shallow subsurface[END_REF][START_REF] Rivière | A set of measures for the systematic classification of the nonlinear elastic behavior of disparate rocks[END_REF]. In certain settings the coseismic deformation may be limited to shearing of microasperities and grain to grain contacts, suggesting that substrate damage and (complete) healing could be described by a model inspired by non-linear elasticity theory and experiments [START_REF] Guyer | Nonlinear Mesoscopic Elasticity: Evidence for a New Class of Materials[END_REF], Ten Cate et al., 2000, Xia and Jonhson, 2005[START_REF] Broda | Modelling of nonlinear crackwave interactions for damage detection based on ultrasound-A review[END_REF], Riviere et al., 2015[START_REF] Snieder | The time dependence of rock healing as a universal relaxation process, a tutorial[END_REF], Sens-Schonfelder et al., 2019). Alternative models, more appropriate for higher strains, could (also) involve processes, such as the opening and closing of cracks and micro-cracks, where crack closing and sealing would constitute only partial recovery. Theory and experiments relating to fracture mechanics may be better suited to explore this in more detail (e.g., Brantut 2015 and references therein).

Last, we note that the type and magnitude of any seasonal dV/V variations are different between stations and also vary at a station for different frequency bands (Fig 3,6). An effort to constrain the links between such pseudo-cyclic velocity variations and the hydrological cycle is an essential step towards a better resolution of the magnitude and evolution of earthquakeinduced velocity changes in the shallow substrate. Such an analysis may refine importantly the best estimate and uncertainty of the Chi-Chi and Gorkha case. Additionally, a well-calibrated hydrological model (Hillers et al., 2014, Clements andDenolle, 2018) may allow using the difference between pre-and post-seismic coupling of dV/V and hydrological measurements to retrieve mechanical changes.

Advancing the field

Standardizing seismological dV/V to overcome the limitations of geomorphic data

Constraints on landslide susceptibility with high temporal resolution are inherently difficult to obtain after earthquakes, because they require frequent, intense rainfall to repeatedly trigger new landslides. Taiwan meets this requirement, but in many regions intense rainstorms are less frequent, and/or topography is less steep, and both factors limit the frequency of landslide events. As a result, in areas with a (semi-)arid climate, such as the epicentral area of the 1994 Northridge, California, earthquake, it may be impossible to discern any seismically-induced change to landslide rates. Even in the Himalaya, where landslides can be frequent during the monsoon season [START_REF] Marc | Long-term erosion of the Nepal Himalayas by bedrock landsliding: the role of monsoons, earthquakes and giant landslides[END_REF], the long, dry period between September and May, means that landslide susceptibility can only be constrained with confidence once per year. This limits our ability to resolve the dynamics of landslide susceptibility recovery, which, based on our limited evidence, appears unlikely to last more than one to a few years [START_REF] Marc | Transient changes of landslide rates after earthquakes[END_REF][START_REF] Marc | Long-term erosion of the Nepal Himalayas by bedrock landsliding: the role of monsoons, earthquakes and giant landslides[END_REF]. For this reason, it is unlikely that we shall ever be able to systematically and precisely constrain the rate at which the landslide susceptibility recovers, with help of geomorphic data. This crucial deficiency underlines the potential utility of high-resolution seismic velocity monitoring if it can indeed be used to quantify and track landscape properties. In such case, estimation of the shallow dV/V coseismic drop at various stations, in the weeks following an earthquake, could allow to determine which zone of the epicentral area, if any, are likely to produce exceptional landsliding in the following rainfalls events. Then early estimates of  S , refined continuously with subsequent recordings, could allow decisionmakers to anticipate when the landslide hazard will have returned to the nearbackground level, and thus to adapt mitigation and reconstruction strategies.

To go beyond temporal co-evolution of landslide susceptibility and seismic velocity, it would be useful to be able to compare the magnitude of both recoverable changes, P L and P S . P L is defined by normalizing for the background landslide susceptibility to rainfall, making it possible to compare areas with different geological and hydrological conditions. A similar normalization has not yet been developed for P S . For P S estimates to scale primarily with subsurface damage, it would be required to (1) have similar uniform seismic interferometry processing, (2) isolate changes in the near-surface only (e.g., 0-100m), and (3) normalize for the medium susceptibility around the station.

Uniform seismic interferometry processing is, in theory, straightforward, but it may be difficult in practice, due to local subsurface noise conditions leading to unstable noise correlation that do not allow efficient retrieval of the seismic velocity at some frequencies, or only with some coda lapse time, as is the case for stations SSLB and GUMBA.

Similarly, isolating near-surface changes requires a multi-frequency analysis, at least up to 2-4Hz, or 4-8Hz and a vertical velocity model. Where boreholes are not available, active or passive seismic experiments may help to constrain the shallow velocity structure, and to adequately quantify the magnitude of shallow changes. New methods focusing on correlating specific types of waves may also allow to better locate shallow changes [START_REF] Mordret | Noise-based ballistic wave passive seismic monitoring -Part 2: surface waves[END_REF].

Previous work has shown that velocity changes are station-dependent [START_REF] Takagi | Coseismic velocity change in and around the focal region of the 2008 Iwate-Miyagi Nairiku earthquake[END_REF][START_REF] Viens | Complex near-surface rheology inferred from the response of greater Tokyo to strong ground motions[END_REF] and that normalization for the station susceptibility, by defining the velocity changes for a reference strain level, could yield a better characterization of disturbances. For this purpose, Earth tides could be used as a universal forcing on the substrate with which station susceptibility can be estimated independently [START_REF] Hillers | In situ observations of velocity changes in response to tidal deformation from analysis of the high-frequency ambient wavefield[END_REF][START_REF] Takano | Seismic velocity changes caused by the Earth tide: Ambient noise correlation analyses of small-array data[END_REF], 2019, Sens-Schönefelder and Eulenfeld, 2019). Addressing these three challenges in future case studies would be an important step toward testing the expected scaling between P L and P S , and toward eventually opening a way for quantitative monitoring of landscape subsurface strength changes with seismic interferometry.

Introducing additional evidence

Before concluding, here we summarize the available constraints on post-earthquake hydro-mechanical landscape changes for different regions of the world, namely the ones analyzed in this study (Japan, Taiwan, Nepal), as well as areas in California and in the Wenchuan epicentral areas (Table 1). We also briefly discuss potential analyses that could be performed to complement some cases.

In addition to an improved standardization of the seismological measurements detailed above, it may also be useful to compare observed seismic and geomorphic changes to hydrological changes, for example using well water level response to tides. Indeed postearthquake increase in permeability has been related to coseismic cracking of the near surface materials [START_REF] Wang | Coseismic release of water from mountains: Evidence from the 1999 (Mw = 7.5) Chi-Chi, Taiwan, earthquake[END_REF][START_REF] Wang | Large earthquakes create vertical permeability by breaching aquitards[END_REF], its magnitude and evolution could relate to landslide susceptibility and dV/V. Illustrating the interest of this approach, the seismological and hydrological recoveries in the epicentral area of the 2008 Wenchuan earthquake have both been estimated to have taken about one year [START_REF] Obermann | Depth sensitivity of seismic coda waves to velocity perturbations in an elastic heterogeneous medium[END_REF], Liu et al., 2014, Xue 2016). In contrast, in Taiwan, a number of wells near the epicenter of the Chi-Chi earthquake also displayed enhanced permeability, but only for 1-2 months [START_REF] Wang | Large earthquakes create vertical permeability by breaching aquitards[END_REF], much shorter than the response we observed for dV/V and landsliding (Fig2), possibly because wells were located in alluvial fans. Understanding the cause of the permeability change and its control remains challenging, but nevertheless, a detailed comparison of the dynamics of such changes with the dV/V or landslide susceptibility may give new insights into coupled mechanisms.

In California, the response of shallow permeability to several earthquakes have been shown to last for months to years and increase with local shaking intensity [START_REF] Elkhoury | Seismic waves increase permeability[END_REF]. Additionally, the nearby seismic velocity response to tidal forcing has been constrained [START_REF] Hillers | In situ observations of velocity changes in response to tidal deformation from analysis of the high-frequency ambient wavefield[END_REF]. In parallel, significant velocity drop and nonlinear relaxation lasting months to years have been evidenced for the San Simeon and Parkfield earthquakes [START_REF] Brenguier | Postseismic Relaxation Along the San Andreas Fault at Parkfield from Continuous Seismological Observations[END_REF], even if the mechanism at play and their depth remain discussed [START_REF] Sleep | Depth of Rock Damage from Strong Seismic Ground Motions near the 2004 Parkfield Mainshock[END_REF][START_REF] Wu | Constraining depth range of S wave velocity decrease after large earthquakes near Parkfield, California[END_REF]. We believe that attempting to couple these analyses, which so far have been performed over different periods or areas (Table 1), could yield robust insights into both the nature of the seismic recovery and the origin of subsurface mechanical changes. However, landsliding is likely insufficiently frequent in most parts of California to allow a direct comparison of seismic and geomorphic responses.

Finally, observational constraints on the recovery of landslide susceptibility may be complemented with observations of other associated processes of mass wasting. This could include the meteorological threshold for the initiation of debris flows (Table 1). A reduction of this threshold has been found to have persisted for about five years after the Chi-Chi earthquake in Taiwan (e.g., [START_REF] Shieh | Variability in rainfall threshold for debris flow after the Chi-Chi earthquake in central Taiwan, China[END_REF] and one year after the Gorkha earthquake in Nepal [START_REF] Dahlquist | Initiation and Runout of Post-Seismic Debris Flows: Insights From the 2015 Gorkha Earthquake[END_REF], consistent with independent estimates of landslide susceptibility changes in these areas [START_REF] Marc | Transient changes of landslide rates after earthquakes[END_REF][START_REF] Marc | Long-term erosion of the Nepal Himalayas by bedrock landsliding: the role of monsoons, earthquakes and giant landslides[END_REF]. Moreover, an increase of the landslide susceptibility lasting up to several years has been inferred from the temporal evolution of the local debris flow threshold in the epicentral area of the 2008 Wenchuan earthquake in China [START_REF] Fan | Earthquake-Induced Chains of Geologic Hazards: Patterns, Mechanisms, and Impacts[END_REF]. However, post-seismic debris flow activity may reflect the evacuation of coseismic landslide deposits rather than any landscape property changes. The relevant temporal evolution of the debris flow threshold should be quantified by separating it from the effect of deposit remobilization (Domenech et al., 2019).

CONCLUSIONS

In this paper we have explored the hypothesis that seismic velocity changes derived from interferometry can be used to track landscape subsurface hydromechanical properties and their response after large earthquakes. These properties are essential to hydrological, geomorphological and natural hazard studies. To test this hypothesis, we have aggregated geomorphological and seismological data that constrain the temporal evolution of the landslide susceptibility and the seismic velocity in the epicentral areas of four intermediate to large, shallow earthquakes (Mw 6.6-7.9). We re-analyzed existing data in Japan, Taiwan and Nepal and added new analyses of velocity changes from ambient seismic noise correlation following the 1999 Chi-Chi and 2015 Gorkha earthquakes in Taiwan and Nepal, respectively.

At present, the Chi-Chi earthquake is the only case in which both phenomena can be shown to co-evolve over about 4 years. For the studied cases in Japan and Nepal the landslide susceptibility and seismic velocity evolution are poorly constrained, but the first-order estimates of their recovery times agree in these three cases. These recovery times were substantially shorter than in Taiwan (a few months to two years). Additional evidence indicates that most of the coseismic change and recovery in Japan may have occurred at very shallow depths of <200 m below the surface. This combined evidence suggests that it may eventually be possible to use seismic interferometry to track the response of near-surface geomaterials in steep landscapes to large earthquakes and to monitor the propensity of such landscapes to slope failure due to subsequent triggers.

We identify several challenges in order to further test our initial hypothesis and to unlock the potential of seismic interferometry to track landscape properties. Efforts are needed to determine velocity changes in several high-frequency bands and to constrain the local subsurface velocity structure to better isolate near-surface changes. Further, and related, efforts are required to normalize the velocity changes relative to depth and material susceptibility in the sensitive area of the seismic stations, so that geomorphic and seismic estimates of the magnitude of the near-surface damage can be compared. Finally, it is important to extend the analysis we have pursued here to other earthquakes for which various geomorphic constraints exist (e.g., the Wenchuan earthquake) and to complement case studies with information on substrate permeability changes. Ultimately, the recovery timescale does not seem to correlate with earthquake magnitude, and its relation to potential drivers, such as the seismic rupture, aftershock sequence or post-seismic deformation, remains to be elucidated. 2016) are classified in three groups: the stations without significant velocity changes at any frequency bands (red), the ones with some variations in one frequency band (blue) and the ones with large changes in at least two frequency bands (yellow). In (C), only two stations where recording before the earthquake, without velocity changes for TDCB, while in SSLB the wavefield coherence only allows the retrieval of velocity changes at high frequencies (2-4Hz). In (E), the epicenter of the main shock and largest aftershock are shown, as well as names of the valley where landslides were mapped. the stations without significant velocity changes at any frequency bands (red), the ones with some variations in one frequency band (blue) and the ones with large changes in at least two frequency bands (yellow). In (C), only two stations where recording before the earthquake, without velocity changes for TDCB, while in SSLB the wavefield coherence only allows the retrieval of velocity changes at high frequencies (2-4Hz). In (E), the epicenter of the main shock and largest aftershock are shown, as well as names of the valley where landslides were mapped. Total landsliding (Area and volume) in Nepal, caused by monsoons in the three valleys affected by the EQ against the landsliding from the non-affected valley (Kali Gandaki, see Fig 1). Numbers correspond to years after 2010, 2015 being the year of the earthquake. When excluding the year 2015, the landsliding seems roughly proportional to 0.7 times the total landsliding in the Kali Gandaki valley, as shown by the solid and dashed lines that are the best linear fit of total volume and total area, respectively. Right: Mean landslide susceptibility through time (Normalized by 0.7 times the Kali Gandaki landslide amount) for landslide area and volume. We only focus on volume as for the other cases.

List of figures captions

Figure 1 :

 1 Figure 1 : (A) Locations of the three study areas. Seismic and GPS station locations relative to

Figure 2 :

 2 Figure 2: Evolution of rainfall-normalized landslide susceptibility (A,) and of the seismic velocity

Figure 3 :

 3 Figure 3: Velocity changes from single-station cross-correlation (red dots) and models for

Figure 4 :

 4 Figure 4: Seismic velocity changes from auto-correlation (Z only) of the signal in the frequency

Figure 5 :

 5 Figure 5: Total landsliding (Area and volume) in Nepal, caused by monsoons in the three valleys

Figure 6 :

 6 Figure 6: Modeled exponential recovery for poorly constrained landslide time-series in the

Figure 7 :

 7 Figure 7: Relaxation and recovery time (estimated as 3 times the relaxation timescales) for the

Figure 1 (

 1 Figure 1 (A) Locations of the three study areas. Seismic and GPS station locations relative to epicenter (star), seismogenic fault plane (dashed polygon with solid line representing the upper end of the fault) and the landslide mapping area (polygon(s)) in the epicentral area of the Niigata 2004 (B), Chi-Chi 1999 (C), Iwate 2008 (D) and Gorkha 2015 (E) earthquakes. In (B) and (D), seismic stations reanalyzed from Hobiger et al. (2016) are classified in three groups: 850the stations without significant velocity changes at any frequency bands (red), the ones with some variations in one frequency band (blue) and the ones with large changes in at least two frequency bands (yellow). In (C), only two stations where recording before the earthquake, without velocity changes for TDCB, while in SSLB the wavefield coherence only allows the retrieval of velocity changes at high frequencies (2-4Hz). In (E), the epicenter of the main shock and largest aftershock are shown, as well as names of the valley where landslides were mapped. 855

Figure 2

 2 Figure 2 Evolution of rainfall-normalized landslide susceptibility (A,) and of the seismic velocity changes (B) relative to a pre-860 earthquake mean. Best-fit of the exponential relaxations are shown in both panels. The grey band in A shows the estimated uncertainty around the background landslide susceptibility. The horizontal dashed black line in B shows the original mean velocity and underlines the partial recovery.

Figure 3

 3 Figure 3Velocity changes from single-station cross-correlation (red dots) and models for various relaxation times (solid lines) for each frequency ranges, 0.5-1 Hz, 1-2 Hz and 2-4 Hz. Considred stations are ICEH (A,C,E) and ICWH (B,D,F), for the 2008 Iwate-Miyagi earthquake, and NGOH (I, K, M) and YNTH (J, L, N) for the 2004 Niigata earthquake (Locations are inFig 1). The ratios of the variance of the best fit model over the variance of other models, used to retrieve confidence intervals with an F-test are in (G), (H), (O) and (P), for either station and each frequency. There, the black curve combines the residuals of all frequency ranges considered in the above plots. The critical variance ratio above which models are indistinguishable from the best model (at a 95% confidence level) is indicated by a black horizontal line.

Figure 4 :

 4 Figure 4:

  Figure 5:

Acknowledgments

This work was initiated during Odin Marc PhD, funded by a fellowship in the EU Marie-Curie Initial Training Network TOPOMOD, project number 264517. Part of this work has been commented on by several anonymous reviewers in previous submissions. The authors have no conflict of interest related to this work.

Data and Ressources

This work also includes three additional figures in an electronic supplemental material. We used seismic waveform data from the Hi-net and BATS networks, provided by the National Research Institute for Earth Science and Disaster Resilience (NIED), Tsukuba, Japan (https://doi.org/10.17598/NIED.0003, last accessed 08/2020), and the Institute of Earth Sciences, Academia Sinica, Taiwan (https://doi.org/10.7914/SN/TW, last accessed 08/2020),

respectively. Laurent Bollinger and the Nepal Seimological Center have provided us with the records from station GUMBA and are warmly thanked. Landslide data can be accessed in the supplementary material of Marc et al., 2019 (For Nepal) or by contacting the corresponding author (for Taiwan and Japan).