Julian Kuehnert
email: kuehnert@ipgp.fr

Anne Mangeney

Yann Capdeville

Jean-Philippe P Métaxian

Luis Fabian Bonilla

Eleonore Stutzmann

E Chaljub

Patrice Boissier

C Brunet

Philippe Kowalski

F Lauret

Clément Hibert

Simulation of Topography Effects on Rockfall-Generated Seismic Signals: Application to Piton de la Fournaise Volcano

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Interactions between seismic wave fields and complex surface geometries can locally modify seismic ground motion. Anomalously strong shaking on hilltops and mountain ridges or flanks, often causing severe structural damage to buildings (W. H. K. [START_REF] Lee | Digital seismograms of selected aftershocks of the Northridge earthquake recorded by a dense seismic array on February 11, 1994 at Cedar Hill Nursery in Tarzana[END_REF][START_REF] Hartzell | Initial investigation of site and topographic effects at Robinwood Ridge, California[END_REF][START_REF] Hough | Localized damage caused by topographic amplification during the 2010 M7.0 Haiti earthquake[END_REF] or triggering earthquake-induced landslides [START_REF] Meunier | Topographic site effects and the location of earthquake induced landslides[END_REF][START_REF] Harp | Relation of landslides triggered by the Kiholo Bay Earthquake to modeled ground motion[END_REF], have been related to seismic amplification due to such topographic effects. Data from field experiments support the assumption of amplified ground motion at the top compared to the bottom of a mountain [START_REF] Davis | Observed effects of topography on ground motion[END_REF]Pedersen et al., 1994;[START_REF] Spudich | Directional topographic site response at Tarzana observed in aftershocks of the 1994 Northridge, California, earthquake: Implications for mainshock motions[END_REF].

Numerous studies have tried to quantify numerically the topographic effect on seismic waves generated by deep sources. [START_REF] Geli | The effect of topography on earthquake ground motion: a review and new results[END_REF] provided an extensive review of previous studies together with new results from more complex models (i.e. including subsurface layering and neighboring ridges). Using an earthquake simulation with three-dimensional topography, [START_REF] Bouchon | Seismic response of a hill: The example of Tarzana, California[END_REF] found that a small hill of less than 20-m high can amplify ground acceleration by 30% to 40% for frequencies between 2 Hz and 15 Hz.

Using the 3D spectral element method, S. J. [START_REF] Lee | Effects of realistic surface topography on seismic ground motion in the Yangminshan region of Taiwan based upon the spectral-element method and LiDAR DTM[END_REF] studied the effects of high-resolution surface topography. They found that values of Peak Ground Acceleration (PGA) can be increased up to 100% relative to simulations on a flat surface and reported an increase in cumulative kinetic energy of up to 200% as a result of increased duration of shaking linked to complex reflection and scattering processes during the interaction of the seismic waves with the topography.

Yet, because of complex patterns of amplification and deamplification, it is difficult to quantify the effect of topography in a generic way. [START_REF] Maufroy | Frequency-scaled curvature as a proxy for topographic site-effect amplification and groundmotion variability[END_REF] proposed to use the topography curvature, smoothed over a characteristic length depending on the studied wavelength, as a proxy for amplification factors. Based on the NGA-West2 earthquake catalog [START_REF] Ancheta | NGA-West2 database[END_REF], [START_REF] Rai | Empirical Terrain-Based Topographic Modification Factors for Use in Ground Motion Prediction[END_REF] showed statistical biases of site residuals in the ground-motion prediction equation (GMPE, [START_REF] Chiou | Update of the Chiou and Youngs NGA model for the average horizontal component of peak ground motion and response spectra[END_REF] towards relative elevation and smoothed curvature and suggested topographic modification factors dependent on signal frequency and relative elevation. In addition to these manuscript submitted to JGR: Solid Earth successful findings, other authors have pointed out the complex coupling between topography and the underlying soil structure that must not be neglected when estimating topographic amplification [START_REF] Assimaki | Ground-motion observations at hotel Montana during the M 7.0 2010 Haiti earthquake: Topography or soil amplification[END_REF][START_REF] Hailemikael | Ground-motion amplification at the Colle di Roio ridge, central Italy: A combined effect of stratigraphy and topography[END_REF][START_REF] Wang | Forward and inverse studies on scattering of Rayleigh wave at surface flaws[END_REF][START_REF] Jeong | How topographydependent are topographic effects? Complementary numerical modeling of centrifuge experiments[END_REF].

All the studies mentioned above investigate topographic effects on a seismic wave field of vertical incidence. S. J. Lee, Komatitsch, et al. (2009) investigates the influence of the source depth on ground motion amplification and demonstrates that amplification in a basin can be reduced when a mountain range is located between the basin and a shallow source. This suggests that surface topography can have a pronounced influence on the propagation of surface waves subjected to an accumulated effect of scattering, diffraction, reflection, and conversion. It is crucial to enhance our understanding of these mechanisms for the study of shallow seismic sources that have gained increasing attention in the emerging field of environmental seismology [START_REF] Larose | Environmental seismology: What can we learn on earth surface processes with ambient noise[END_REF]. Several authors have investigated numerically the interaction of surface waves with 2D surface geometries such as corners, hills or canyons [START_REF] Munasinghe | Finite difference analysis of rayleigh wave scattering at vertical discontinuities[END_REF][START_REF] Weaver | on Diffuse Waves in Solid Media[END_REF][START_REF] Sánchez-Sesma | The influence of topography on the propagation and scattering of surface waves[END_REF][START_REF] Sánchez-Sesma | The influence of topography on the propagation and scattering of surface waves[END_REF]Zhang et al., 2018;[START_REF] Wang | Forward and inverse studies on scattering of Rayleigh wave at surface flaws[END_REF]. [START_REF] Ma | Effects of large-scale surface topography on ground motions, as demonstrated by a study of the San Gabriel Mountains, Los Angeles, California[END_REF] demonstrated that a topographic feature 10 times smaller than the wavelength can still considerably reduce the amplitude of by-passing surface waves. Similar to S. J. Lee, Komatitsch, et al. (2009), they simulated the shielding effects of large-scale topography on fault-generated surface waves using a 3D model of the San Gabriel Mountains, Los Angeles, California, finding amplification factors in peak ground velocity (PGV) of up to +50% on the source-side of the mountain range and up to -50% on the opposite site. L. [START_REF] Wang | Effect of near-surface topography on high-frequency Rayleigh-wave propagation[END_REF] modeled the influence of an uplifted and a depressed topography on the wave field. Comparing amplitudes and frequency content between source side and far-source side, they found that the depressed topography caused stronger contrasts than the uplifted topography, especially for steeper slopes and at higher frequencies.

The present study is focused on seismic waves generated by rockfalls. Different than the source mechanism of earthquakes, rockfall seismic sources can generally be described by impulse forces on the Earth's surface. Seismic signals from rockfalls, or more generally from landslides, have been demonstrated to be very useful to classify and locate events as well as constrain flow dynamics and rheology (e.g. [START_REF] Vilajosana | Rockfall induced seismic signals: Case study in Montserrat, Catalonia[END_REF][START_REF] Deparis | Analysis of rock-fall and rock-fall avalanche seismograms in the French Alps[END_REF][START_REF] Favreau | Numerical modeling of landquakes[END_REF][START_REF] Hibert | Slope instabilities in Dolomieu crater, Réunion Island: From seismic signals to rockfall characteristics[END_REF][START_REF] Dammeier | Characterization of alpine rockslides using statistical analysis of seismic signals[END_REF][START_REF] Moretti | Numerical modeling of the Mount Steller landslide flow history and of the generated long period seismic waves[END_REF][START_REF] Bottelin | Seismic and mechanical studies of the artificially triggered rockfall at Mount Néron (French Alps, December 2011)[END_REF]. However, as landslides predominantly occur in areas of strong topographic relief, the measurements can be strongly influenced by topography variations leading to erroneous landslide estimates. For example, to calculate landslide volumes, the generated seismic energy is estimated from seismic recordings [START_REF] Hibert | Slope instabilities in Dolomieu crater, Réunion Island: From seismic signals to rockfall characteristics[END_REF]. At the same time, energy estimations can vary from station to station. The present work shows that the topography studied here can partly explain amplitude variations between seismic stations.

In the following, after introducing the study site located at Dolomieu crater on Piton de la Fournaise volcano, Reunion Island, the numerical model for the SEM simulations is defined, entailing a discussion on the seismic velocity profile at Piton de la Fournaise.

As the mesh size affects the computational cost, different topography resolutions are compared. Then, topography induced amplification is computed, depending on the underlying velocity model, by means of peak ground velocity (PGV) and total kinetic energy for both vertical and horizontal sources. In an attempt to quantify the dependencies on geometric parameters, the influence of variations in crater depth and curvature are investigated.

Finally, real seismic signals generated by rockfalls at Dolomieu crater are analyzed.

Simulated and observed inter-station spectral ratios are compared, making it possible to examine the spectral content of the signals independently of the rockfall source. Additionally, the seismic signature of single rockfall impacts is investigated. To compare signal amplitude and frequency content between observations and simulations, impact force and duration are estimated based on Hertz contact theory [START_REF] Hertz | Die lieber die Berührung fester elastischer Körper[END_REF].

The study site is located on Piton de la Fournaise volcano, Reunion Island, presented in Figures 1a and1b. Its summit is characterized by the 340 m deep Dolomieu crater that collapsed in 2007 (e.g. [START_REF] Staudacher | The April 2007 eruption and the Dolomieu crater collapse, two major events at Piton de la Fournaise (La Réunion Island, Indian Ocean)[END_REF]. Because of the instability of the crater walls, rockfall events are frequently observed within the crater [START_REF] Hibert | Slope instabilities in Dolomieu crater, Réunion Island: From seismic signals to rockfall characteristics[END_REF][START_REF] Hibert | Automated identification, location, and volume estimation of rockfalls at Piton de la Fournaise volcano[END_REF], 2017;[START_REF] Durand | On the Link Between External Forcings and Slope Instabilities in the Piton de la Fournaise Summit Crater, Reunion Island[END_REF][START_REF] Derrien | Multi-temporal airborne structure-from-motion on caldera rim: Hazard, visitor exposure and origins of instabilities at Piton de la Fournaise[END_REF].

The high quantity of events together with a dense seismic network monitored by the Observatoire Volcanologique du Piton de La Fournaise (OVPF) provide excellent conditions for the study of rockfalls. Using recorded seismic signals, past studies have investigated the links between rockfall activity and external forcings such as rain or seismicity, the spatio-temporal evolution of rockfall occurrences as well as their volumes [START_REF] Hibert | Dynamics of the Bingham Canyon Mine landslides from seismic signal analysis[END_REF][START_REF] Hibert | Spatio-temporal evolution of rockfall activity from 2007 to 2011 at the Piton de la Fournaise volcano inferred from seismic data[END_REF][START_REF] Durand | On the Link Between External Forcings and Slope Instabilities in the Piton de la Fournaise Summit Crater, Reunion Island[END_REF].In addition to the seismic stations, three cameras positioned on the crater rim continuously monitor rockfall activity. This makes it possible to correlate video images and rockfall seismic signals.

For example, Figures 1c-1e show images and seismic signals of a rockfall on the southern crater wall on February 28, 2016. The rockfall consisted of mainly three boulders that can clearly be traced on the video provided in the Supporting Information. Each boulder took around 30 s to move from the top to the bottom of the crater.

The first movement can be seen in snapshot i). At that time, a large signal amplitude was recorded on station DSO, located very close to the source position. Subsequently, the rockfall traveled through a small valley in ii) and accelerated towards the position in iii). The acceleration of the boulder resulted in strong impacts which can be seen on both the signal and the spectrogram after time iii) at all stations. At the time corresponding to snapshot iv), the first boulder arrived at the crater bottom, whereas a second boulder was half-way down. Again, strong amplitudes were measured around time iv), probably corresponding to the second boulder. Around time v), the last movements of a third block are visible. After this, residual granular activity distributed on the flank can be observed on the video. Signal amplitudes decay accordingly.

Note that station DSO recorded very strong signals in the beginning, while signal amplitudes increased slowly at the other stations. This is certainly related to the changing source-receiver distance. Additionally, topography may have influenced the signal amplitudes depending on the source position relative to the receiver position. From the spectrograms, we can see that the main frequency content was between 3 Hz and 20 Hz.

SEM simulations

To study the effect of topography on rockfall seismic signals recorded at different stations, the seismic wave propagation was simulated using the 3D Spectral Element Method (SEM, e.g. [START_REF] Festa | The Newmark scheme as velocity-stress time-staggering: An efficient PML implementation for spectral element simulations of elastodynamics[END_REF][START_REF] Chaljub | Spectral-element analysis in seismology[END_REF]. The seismic impulse response was modeled by implementing a point force at the surface of the domain in the form of a Ricker wavelet with a dominant frequency of 7 Hz, covering a bandwidth from 2 Hz to 20 Hz, which is predominantly observed for rockfalls at Dolomieu crater. The source magnitude was set to unity except for analysis of a single impact rockfall in section 5.3.

Moving rockfall source positions and poorly known subsurface properties require many simulations with different configurations. For easy reference, a table is provided in Appendix A, listing the simulation configurations used in the different sections of the article. absorbing PML boundaries (Perfectly Matched Layers, e.g. [START_REF] Festa | The Newmark scheme as velocity-stress time-staggering: An efficient PML implementation for spectral element simulations of elastodynamics[END_REF] are added on the sides and bottom of the domain. The elements are successively deformed in the vertical direction to accommodate the topography provided by a Digital Elevation Model (DEM) with a 10 m resolution. To decrease computational costs, the element size is increased from 10 m to 30 m at 150 m below the surface (Zone of refinement), resulting in a total of 915,704 elements. To filter out short wavelength variations of the fine mesh that cannot be represented in the coarse mesh, a smooth Buffer layer, provided by a low-pass filtered topography, is used as an additional boundary 100 m below the surface. Simulations are implemented using a polynomial degree of 5, i.e. 6 GLL points (Gauss-Lobatto-Legendre, e.g. [START_REF] Chaljub | Spectral-element analysis in seismology[END_REF] per element in each direction.

Mesh of the Earth model

Given the high computational costs (i.e. CPU times of up to 460 days, with 10 cores per CPU, for the heterogeneous velocity model), a mesh with a reduced topography resolution of 20 m is used in the first part of this study in which different velocity models are explored with a fixed point source located at the southern crater wall. This mesh is built using elements with a constant side length of 20 m, which reduces the total number of elements to 550,000, accordingly decreasing the CPU time to 145 days for the heterogeneous velocity model. The reduced number of elements also increases memory efficiency when displaying snapshots.

In the second part of the study, when comparing simulations to observations of rockfalls at Dolomieu crater, the model with high-resolution topography is used. Here we apply the reciprocity principle (Bettuzzi, 2009), i.e. the synthetic source is located at the position of the real seismometers (BON, BOR, DSO and SNE) and the wave field is recorded on a 10×10 m grid of stations across Dolomieu crater. In this way, the impulse responses of all potential rockfall sources are modeled with just one simulation per seismometer and per channel.

The mesh of the flat-surface reference model is built with elements with 20 m side lengths. Cross-sections through all meshes (i.e. those with a flat surface; with 20 m-resolution topography; with 10 m-resolution topography) are provided in the Supporting Information.

Velocity model

Three different velocity models are used: (1) a homogeneous model, (2) a model with shallow low S-wave velocity layer, and (3) a model with smoothly increasing velocity as proposed by [START_REF] Lesage | A generic model for the shallow velocity structure of volcanoes[END_REF] for shallow volcano structures. The velocity-depth profiles are illustrated in Figure 2b and summarized in Table 1. The generic model by [START_REF] Lesage | A generic model for the shallow velocity structure of volcanoes[END_REF] is based upon measurements on multiple andesitic and basaltic volcanoes. P-and S-wave speeds c i are expressed as follows:

c i (z) = c i0 [(z + a i ) αi -a αi i + 1], i = P, S, (1)
where z is the depth below the surface, c i0 are the velocities at zero depth, and α i and a i are fitting parameters as defined in Table 1.

The velocity profiles are compared to the S-wave velocity model inverted from ambient noise recordings at Piton de la Fournaise by [START_REF] Mordret | Three-dimensional shear velocity anisotropic model of Piton de la Fournaise Volcano (La Réunion Island) from ambient seismic noise[END_REF]. The orangeshaded zone shown in Figure 2b corresponds to depth-profiles extracted from the inverted 3D model in the vicinity of Dolomieu crater. Good agreement is observed with the Lesage generic velocity profile. The discrepancy in the first 100 m can be caused by missing highfrequency content in the model of [START_REF] Mordret | Three-dimensional shear velocity anisotropic model of Piton de la Fournaise Volcano (La Réunion Island) from ambient seismic noise[END_REF], who inverted frequencies below 2.5 Hz.

In order to further validate the Lesage generic model for our study site, Rayleigh velocity dispersion curves from noise measurements at a mini-array located around station BON are compared in Figure 2c with theoretical dispersion curves of the Lesage generic model. Picks from the mini-array measurements are determined using the Modified Spatial Autocorrelation (MSPAC) Toolbox [START_REF] Köhler | Assessing the reliability of the modified three-component spatial autocorrelation technique[END_REF][START_REF] Wathelet | Array performances for ambient vibrations on a shallow structure and consequences over Vs inversion[END_REF] 2a. There are two options when implementing a velocity-depth profile on a 3D numerical domain with topography. The first possibility is to keep the velocity laterally homogeneous and excavate a surface corresponding to the topography. The second possibility is to adjust the velocity profile vertically so that it follows the topography elevation. Either way, the subsurface velocity structure is influenced, unless it is homogeneous. We chose the second option that we believe is geologically more reasonable because a main cause of velocity variation is the compaction of material with depth due to increasing overburden pressure.

Rock density ρ as well as quality factors Q P and Q S for intrinsic attenuation of P-and S-wave velocity, respectively, are chosen based on previous studies on Piton de la Fournaise and similar volcanoes [START_REF] Battaglia | Location of seismic events and eruptive fissures on the Piton de la Fournaise volcano using seismic amplitudes[END_REF][START_REF] O'brien | Volcano topography, structure and intrinsic attenuation: Their relative influences on a simulated 3D visco-elastic wave-content[END_REF][START_REF] Hibert | Slope instabilities in Dolomieu crater, Réunion Island: From seismic signals to rockfall characteristics[END_REF]. All parameters are summarized in Table 1.

This work focuses on the topography effect, but it is important to have an idea of the effect of 3D-medium heterogeneities. Difficulties arise in this respect as there is lack of knowledge on the distribution of heterogeneities that is hard to invert from seismograms alone [START_REF] Imperatori | Broad-band near-field ground motion simulations in 3-dimensional scattering media[END_REF]. Nevertheless, a first attempt to simulate scattering effects is made by adding a spatially random velocity perturbation to the Lesage generic velocity model. The magnitude of the velocity deviation reaches 43% and is defined by a normal distribution with a standard deviation of 10%, see Supporting Information.

Topography resolution

The influence of topography resolution on the simulated wave field is investigated to study influences from sub-wavelength topography variations and assess the trade-off between increased resolution and computation costs. Synthetic seismograms of the vertical component are compared (Figure 3a), obtained from models with a flat surface and topography resolutions of 20 m and 10 m. A vertical point force in the form of 7 Hz Ricker source-time function is placed on the southern crater wall, corresponding to rockfall starting position R1 in Figure 1b.

The single-impact source produces a long wave-train of body waves and multiplemode Rayleigh waves. The maximum amplitudes decrease for the models with topography compared to the simulations with the flat model. The large-scale crater topography can redirect the wave-field and cause shadow-zones, as shown in section 4 where the spatial distribution of amplification or deamplification along the surface will be analyzed.

Furthermore, topography causes prolonged and more complex waveforms. For the flat manuscript submitted to JGR: Solid Earth This assumption is supported by the spectra recorded at station BON, Figure 3b.

Differences between the two models with topography become evident above roughly 5 Hz.

This corresponds to a minimum wavelength of 116 m for the fundamental Rayleigh wave (λ ≈ 580 m s -1 ÷5 Hz ≈ 116 m). If we conclude that wavelengths below 116 m are still sensitive to the change in topography resolution, then 1 st -mode Rayleigh waves of above 7 Hz are affected (λ ≈ 800 m s -1 ÷7 Hz ≈ 114 m). This analysis suggests surface waves are sensitive to changes in topography resolution that are 5 times smaller than their wavelength.

The decrease in amplitude at all stations for the high-resolution topography (10 m)

relative to the lower-resolution topography (20 m) suggests that more energy is scattered during the propagation of the surface waves and possibly lost within the subsurface. Interestingly, S. J. [START_REF] Lee | Effects of realistic surface topography on seismic ground motion in the Yangminshan region of Taiwan based upon the spectral-element method and LiDAR DTM[END_REF] found the opposite when comparing waveforms on different topography resolutions for a source deep below the surface. This implies that the source position plays a major role in the effect of topography. On one hand, topography can increase ground shaking and thus trap energy close to the surface. On the other hand, in the case of waves traveling along the surface, the topography can increase scattering and thus prevent energy propagation. Similar conclusions were drawn by S. J. Lee, Komatitsch, et al. (2009), who investigated how topography effects are modulated by the source depth in regard to ground motion in a basin located behind a mountain range.

Wave propagation from a vertical surface load

To better understand wave propagation along the topography and the influence of the subsurface geology, snapshots of the wave field were examined. Simulations for all three velocity models were carried out on the domain with 20 m Dolomieu topography resolution and a vertical point force on the southern crater wall. For the simulation with the homogeneous domain (left column of Fig. 4), we can identify in the first snapshot at time t = 0.8 s the P-wave traveling downwards as being the fastest wave with propagation direction parallel to the shown vertical ground velocity. At time t = 1.6 s, the original S-wave is visible on the bottom of the cross-section.

The S-wave can be identified because the direction of propagation is perpendicular to the vertical ground velocity. Just above, note the newly created S-wave (annotated as RS) that separated at the bottom of the crater from the Rayleigh wave because of the convex topography. Yet, part of the energy continues as a Rayleigh wave along the topography towards the rim of the crater. Also visible is a diffracted surface wave (annotated as Rd). It split from a wave front traveling towards station BOR and took a curved path along the flank of the crater. At time t = 2.0 s we can see this diffracted Rayleigh wave continuing outside the crater and arriving at station BON at a different azimuth than that of the Rayleigh wave that traveled diagonally across the crater and its rim (annotated as Rf). The energy of Rayleigh wave Rf was partly reflected at the crater rim so that a new Rayleigh wave Rr traveled backwards through the crater. Up front (on the far right of the domain), a direct S-wave hits the surface and is partly reflected and converted to build a straight P-wave front traveling downwards at an oblique angle to the horizontal (annotated as SP).

Adding a low S-wave velocity layer (middle column in Fig. 4) drastically changes the wave field because of reflections within this layer and the dispersive character of Rayleigh A1. The absence of intrinsic attenuation for enhanced visibility of wave propagation caused reflections from the boundary on the left at later times.

waves. Looking at the first 2.5 s of the synthetic seismograms, we observe a wave-train with a dispersive character overlaid by multiples. Compared to the homogeneous model, it is more complex and has a longer duration because of internal reflections within the low-velocity layer. At around t = 2.6 s, the waves hit the crater rim opposite the source and are partly reflected just like in the homogeneous case. The snapshots at times t = 2.6 s and t = 3.8 s show, in contrast to the homogeneous case, a much more scattered wave field of irregular amplitude patterns. Similar to S. J. [START_REF] Lee | Effects of realistic surface topography on seismic ground motion in the Yangminshan region of Taiwan based upon the spectral-element method and LiDAR DTM[END_REF], who found characteristic patterns dependent on the resolution of the imposed topography, the characteristic length of these patterns is likely to be related to the resolution of the topography and the flat element surfaces with 20 m side lengths.

For the Lesage generic velocity model (right column of Fig. 4), the majority of energy stays close to the surface of the domain because of the velocity gradient. Scatter-manuscript submitted to JGR: Solid Earth ing of the wave field over the topography is even greater than for the case with a lowvelocity layer and the duration of shaking is increased. From synthetic seismograms (top right of Fig. 4), we can still identify the outward propagation of energy as well as the reflection of part of the energy at the crater rim. The analysis of the simulations shows that a single impact can produce a complex wave field caused by the surface topography and the underlying velocity model.

Regarding scattering from 3D medium heterogeneities, results from simulations in which the Lesage velocity model is randomly perturbed by heterogeneities with a standard deviation of 10% (with maximum excess of 43%) show barely affected synthetic signals, see Supporting Information, indicating that the effective medium is not changed significantly. This example does not prove that scattering is weak as it is possible to design a distribution of characteristic correlation lengths that completely changes the waveforms. However, studies at Dolomieu crater lead us to believe that this is not the case; e.g. [START_REF] Hibert | Slope instabilities in Dolomieu crater, Réunion Island: From seismic signals to rockfall characteristics[END_REF] found that rockfall seismic signals do not exhibit a coda but that their duration corresponds to the rockfall propagation time on videos and Kuehnert et al. ( 2020) tracked rockfall trajectories using simulated inter-station energy ratios with the smooth Lesage velocity model.

Influence of topography on simulated wave propagation

Seismic amplitudes carry crucial information on the seismic source and can be used to infer source locations and acting forces. However, as can be concluded from the simulated wave propagation above, topography together with the underlying geology can strongly influence ground motion. Consequently, the measured amplitudes have to be interpreted according to both source properties (including the resulting radiation patterns) and propagation effects. In the following, topography induced amplification is quantified for different velocity models and different source directions. This can be helpful to better interpret the spatial distribution of amplitudes and eventually account for amplified signals.

Amplification for a vertical source

In order to quantify topographic ground motion amplification, simulations on a model with topography are compared to a reference model with a flat surface. The comparison is performed for both vertical peak ground velocity PGV z and total kinetic energy E. Quantifying PGV amplification is important when interpreting signal amplitudes. However, it does not measure the increased complexity and duration of recorded waveforms caused by scattering and diffraction of the wave field along the topography. These effects can be incorporated by calculating energy amplification. Also, frequency dependencies are not considered. For this reason we will later look at different frequency bands or determine spectral ratios when analyzing observed rockfall signals.

To quantify vertical PGV amplification, the maximum vertical ground velocity is measured at each point on the surface defined on a grid with 30 m spacing. The top row of Figure 5 shows the peak ground velocity ratios PGV z,T /PGV z,F between the three velocity models with topography and the flat reference model.

Similarly, energy amplification is calculated at each grid point by the ratio E T /E F between the models with topography and the flat reference model, where E i (with i = T, F ) is a proxy of the total kinetic-energy density, defined as the square of the recorded ground velocity v, integrated over the total signal duration d:

E i = d v 2 x,i (t) + v 2 y,i (t) + v 2 z,i (t) dt, (2)
The resulting energy amplification is shown in the bottom row of Figure 5 for the three different velocity models. 7, 8 in Table A1. The yellow star illustrates the source position and green triangles mark station locations. Annotations give ratios measured at the station locations as well as the percentage of topographic amplification. Neighboring contour lines differ by 60 m in elevation.

PGV amplification

Analyzing the PGV amplification shown in the top row of Figure 5, the homogeneous model shows a contrast between source side of the crater and the opposite side:

PGV is amplified on the source side and strongly deamplified on the far side. The amplification on the source side (+12% at DSO) can be explained by the simultaneous arrival of surface and direct waves emitted from the source. Deamplification on the farside of the source (-83% at BON and -87% at SNE) can be interpreted as a shadow zone behind the crater related to the diversion of a major part of wave energy downwards into the subsurface because of the crater shape.

For the model with the low-velocity layer, general amplification on the source side and deamplification on the far-source side of the crater are still present but less pronounced (deamplification at station SNE is reduced to -67%) and patterns become more complex (DSO is now deamplified by -19%). The introduction of a low-velocity layer causes more energy to stay at the surface and thus reduces the shadow zone behind the crater.

The uneven topography together with the underlying low-velocity layer causes complicated reflections and wave conversions which lead to increased complexity of amplification patterns.

The contrast between source side and far-source side of the crater decreases further for the Lesage generic velocity model (-45% at DSO, -62% at BON and -35% at SNE). As can be seen on the wave propagation snapshots in Figure 4, the gradient causes energy to stay close to the surface. Whereas a lot of energy is lost downwards because of the crater topography in the homogeneous model as well as in the low-velocity layer model, the velocity gradient in the Lesage generic model guides waves efficiently along the crater topography or back to the surface, which causes a more homogeneous amplification pattern. Scattering away from the surface due to surface roughness as well as conversion from vertical to horizontal energy leads to an overall deamplification in vertical PGV. Still, because of focusing mechanisms of the 3D topography, ray-shaped zones of PGV amplification originating at the source can be observed.

Energy amplification

In general, the amplification patterns of kinetic energy (bottom row of Figure 5)

show more contrast than the PGV ratios. This is because topography does not only influence peak amplitude, but also the complexity and duration of the signal. For the homogeneous model, amplification increases to +41% at DSO and decreases to -92% at BON. Behavior for the model with the low-velocity layer is very similar. For the Lesage generic model, the ray-shaped zones of amplification are considerably more pronounced than for the case of PGV amplification. Given that horizontal ground velocity is considered when computing the kinetic energy, this observation suggests that topography guides both vertical and horizontal energy along the same paths. Note also the increased amplification at parts of the crater cliff ridge, possibly due to the discussed reflection of Rayleigh waves at these positions.

Changing the velocity model modifies the wavelengths, which presents an alternative explanation for the observed differences in the amplification patterns. This explanation was discarded after verifying that the differences still remain for bandpass filtered results, comparing amplification patterns from the different velocity models for coinciding wavelengths as done in Appendix B.

Amplification for a horizontal source

Only vertical surface loads are considered above. However, the rockfall generated basal forces on the ground can also have horizontal components. Here we show amplification patterns for a horizontal source using the Lesage generic velocity model. The directionality patterns are no longer visible when analyzing the amplification of total kinetic energy. This is because all components of the measured ground velocity are considered in the energy calculation. The energy amplification pattern is similar to the one for the vertical source as shown above in Figure 5. This suggests that topography guides seismic energy on trajectories along the surface that mainly depend on the source position rather than the source direction. We will further discuss this hypothesis later when studying inter-station spectral ratios of real rockfall signals.

Surface roughness and crater geometry

The amplification patterns observed in the previous section are characterized by complex spatial distributions. We will now perform tests on domains with synthesized surface topographies in order to better understand the contributions of certain geometric features to the amplification pattern. More concretely, we will study a planar surface with natural roughness as well as synthetic crater shapes of different depths and curvatures. Surface roughness and crater dimensions are defined to resemble our study site on Piton de la Fournaise volcano. The initial domain is a cube of size 2360 m×2360 m× We now investigate the effect of surface topography at different frequency bands.

For this, we quantify as before the amplification of total kinetic energy with respect to the flat reference model. Note that here we present energy instead of PGV as it accounts for both amplitudes and prolonged ground shaking and hence gives a more general picture. Figure 7b shows energy amplification on both the rough planar domain and the domain with a synthetic crater in the 3-7 Hz and 13-17 Hz frequency bands.

We can see that both these frequency bands are influenced by the rough planar surface. As already indicated, the rough topography is band-pass filtered at corner wavelengths 40 m and 100 m and the minimum and maximum wavelengths of fundamental Rayleigh waves are below and above the range of topography wavelengths, respectively. A1.

We observe ray-shaped zones of amplification which are blurred in the lower-frequency band and become sharper towards higher frequencies, because of the shorter interfering wavelengths. The variations of topography seem to guide energy along these ray paths.

In contrast, some areas of pronounced topography variation (visible by the densification of contour lines) seem to shield the propagation of energy and cause shadow zones behind them. This can for example be observed in the north-east direction of the source.

Analyzing the energy amplification on the domain with the synthetic crater, we recognize similarities with the amplification patterns of the previously analyzed planar rough surface. This is because the same surface roughness is used and its imprint is now superimposed on the amplification caused by the crater topography. Globally, the wave field is deamplified behind the crater (as seen from the source position). Higher frequencies seem to be more affected by this than lower frequencies. Nonetheless, even at high frequencies, paths of amplified energy can traverse the crater. This phenomenon seems to be caused by a coupled effect from small-scale (i.e. roughness) and large-scale (i.e. crater) topography variations and is similar to the amplified ray-paths observed on the Dolomieu crater topography (compare to Figure 5). A1.

The energy ratios from the simulations on the domains with synthetic crater shapes are shown on the right-hand side of Figure 8 for the whole frequency range. The amplification pattern varies slightly going from small depth to big depth (Figure 8a). The biggest change is observed behind the crater directly opposite the source. Amplification decreases manuscript submitted to JGR: Solid Earth at this point with increasing crater depth. In contrast, increased amplification is observed inside the crater. These changes in the amplification patterns might be related to interference caused by the symmetric crater form. Going from weak to strong curvature (Figure 8b), the shadow zone behind the crater increases more strongly. This is not only true directly opposite the source position but also diagonally across the crater, suggesting that the increased crater curvature shields off more energy by reflecting or deflecting the wave field sideways or downwards into the subsurface.

The analyses suggest that variations of curvature have stronger effects on ground motion than variations of crater depth for models tested. It is important to note that the wave field is influenced by topography features of scales both below and above the seismic wavelength. This was observed for the experiments with a planar rough surface as well as with the synthetic crater with dimensions (∼ 800 m diameter, ∼ 300 m depth)

largely exceeding the seismic wavelengths.

The experiments on the synthetic model surfaces explored effects of individual aspects of the topography on the wave field. The insights acquired can be transferred to our study side at Piton de la Fournaise volcano as similar scales were chosen deliberately.

Having said this, note that the overall effect on the wave field is governed by the whole configuration and cannot be reduced to an individual feature (e.g. only small-scale versus only big-scale topographic variations). At the same time, the relative position of the source and receiver plays a defining role, i.e. topography related amplitude modifications at a given station can only be predicted if the source position is known.

Seismic signals from rockfalls at Dolomieu crater

We will now study observed seismic signals generated by rockfalls at Dolomieu crater.

As the influence of the topography changes with the source position, we analyze the signals at specific times corresponding to specific rockfall positions. First we will investigate spectral ratios between stations of time windowed rockfall signals. The objective is to clarify as to whether simulations can reproduce the observed spectral ratios when taking into account topography. Subsequently we will focus on a single block impact, identifying its seismic signature and comparing observed and simulated amplitudes by estimating the generated impact force using Hertz contact theory.

Observed spectral ratios between stations

For this analysis, we select three rockfalls with similar trajectories on the southern crater wall corresponding to rockfall location 1 in Figure 1b. The trajectories of the rockfalls were identified from camera recordings. Snapshots of the three events are shown in Figure 9 together with an image of the whole trajectory reconstructed from differences of successive snapshots. Below, the corresponding seismic signals recorded at stations BON, BOR, DSO and SNE are presented.

Station DSO shows the strongest amplitudes, especially in the beginning of the rockfall. This is because the three rockfalls start very close to this station. BON contains the smallest amplitudes, being the furthest station and on the opposite side of the crater.

The dynamics of the three events are not entirely identical. Event 1 consists of a single boulder bouncing down towards the bottom of the crater while other blocks follow with a time lag of around 15 s. In contrast, event 2 consists of two blocks closely following each other down with a time lag of only 4 s, as can be seen on snapshot 2b. Event 3 consists of a main boulder with a smaller block following much later with a lag of about 50 s.

Despite these differences, we compare spectral ratios between stations in time windows R1, R2 and R3 during which the main blocks moved within identical areas. The spectral ratios are computed from the measurements at stations BOR, DSO (vertical component only) and SNE with respect to station BON (note that BON is selected as the reference station as it turns out it is the least affected by local site effects). In order to avoid spurious fluctuations, the spectra are smoothed as proposed by Konno and Ohmachi components. Note the similar behavior of the spectral ratios for each of the events and for each component.

Comparison to spectral ratios from noise recordings (TW-N), from the beginning of event 1 (TW-C1), and from a rockfall that occurred at a different position in the crater (TW-C2, corresponding to trajectory 2 in Figure 1b), shows partly strong deviations from manuscript submitted to JGR: Solid Earth curves R1, R2 and R3. This provides evidence that the spectral ratios are indeed characteristic of the position of the rockfall seismic source. The same analysis is carried out in Appendix D for rockfalls in the southwest, leading to the same conclusion.

Comparison of observed and simulated spectral ratios

The seismic source of a rockfall can be very complex as multiple impacts of different magnitude can occur simultaneously at different positions. Hence, it is very difficult to correctly simulate the rockfall seismic signal, especially at high frequencies. For this reason, spectral ratios between stations are very convenient to compare real and synthetic signals. In this way, the signature of the source is removed from the signal and solely propagation path effects remain. Nevertheless, we have to keep in mind two points when comparing observed and simulated spectral ratios.

Firstly, local subsurface heterogeneities can modify recorded amplitudes and thus influence inter-station ratios. These geological site effects are not considered in the sim- Secondly, different source directions cause different radiation patterns. This is illustrated in Figure 11a, where a force on a flat surface is polarized in vertical (top) and horizontal (bottom) directions. If the radiation pattern is not radially symmetric, which is only the case for vertical ground motion from a vertical source, the spectral ratios are affected depending on the azimuthal position of the respective receivers. The direction of a rockfall seismic source depends both on the rockfall dynamics and on the underlying slope. The generated forces from a boulder impact are schematically illustrated in Figure 11b. The resulting force F r is composed of a force F n normal to the slope and a force F t tangential to the slope, which depends on the slope angle, the direction of movement and the friction between the moving mass and the ground.

In order to analyze the influence of the source direction on the spectral ratios, we compare a vertical force to a normal force and a tangential force. Note that we assume that the tangential force is parallel to the slope of steepest descent. To consider a spatially distributed source in the simulations, the mean spectral ratio is calculated from a selection of multiple sources. This makes it possible to simultaneously evaluate the sensitivity of the curves to the source positions. Seven source positions are picked from a grid with 10 m spacing (see Fig. 1b, picked sources). The area corresponds to the region in which rockfalls 1, 2, and 3 are present during time windows R1, R2, and R3, respectively (see Fig. 9). As the slope dips northwards, the tangential force is orientated in the northdirection. Station BOR is located west of the source position, which is transverse to the source direction. For this reason, a smaller signal amplitude is measured at station BOR in comparison with station BON (ratio < 1), even though BOR is slightly closer to the source. Nevertheless, the tangential force also contains a vertical component that ensures that the ratio is of the same magnitude as the observed ratios. This is different for the spectral ratios of the north component (middle), where a striking discrepancy of more than one order of magnitude results from the vertical source. A vertical force does not generate horizontal transverse energy which is why almost no signal is recorded on the north component at station BOR located eastwards. For the east spectral ratios (right), the tangential force again shows the strongest deviation for reasons similar to those for the vertical component spectral ratios.

Simulated spectral ratios for a model with a flat surface

Simulated spectral ratios for a model with topography

As opposed to the model with a flat surface, a model with the Dolomieu crater topography results in a good agreement between simulated and observed spectral ratios (see Figure 12). Furthermore, very similar values can be observed when comparing the simulations with different source directions, especially towards higher frequencies. This indicates that the spectral ratios are in this case not dominated by the direction of the A1.

source (and the corresponding produced radiation pattern) but rather by propagation along the topography.

Greater deviations between the different source directions are found at lower frequencies, such as for example on the north component of ratio BOR/BON below 3 Hz.

Assuming fundamental Rayleigh waves, this corresponds to wavelengths above 250 m (λ ≈ 750 m s -1 ÷ 3 Hz). With a distance of around 500 m between the source position and station BOR, it is likely that these low-frequency waves have not traveled enough wavelengths in order to be completely dominated by propagation along the topography.

Analyzing the sensitivity of the ratios to the source position, generally larger standard deviations (shaded zone of uncertainty around the mean) are present after intro-manuscript submitted to JGR: Solid Earth A1.

ducing topography compared to the results for the flat model in Figure 11c. This means that a slight change of source position allows more variability of the ratios when considering topography and can eventually better explain the observed spectral ratios.

Clearly, the spectral ratios also depend on the relative source-receiver distance. For example, the high values of ratio DSO/BON result from the fact that the source is very close to station DSO. Furthermore, the values increase towards higher frequencies. This is related to the attenuating properties of the medium that cause the amplitudes of higher frequencies to decrease faster with the distance traveled than lower frequencies.

The analysis suggests that the spectral ratios are characteristic of the source position and dominated by propagation along the topography rather than by the radiation patterns caused by the source directions. To further validate this hypothesis, the same comparison between observations and simulations is carried out in Appendix D for rockfalls located on the southwestern crater wall.

manuscript submitted to JGR: Solid Earth

Seismic signature of a rockfall impact

We will now analyze in detail the seismic signal generated by the single impacts of a rockfall at Dolomieu crater. Interpretation of the signal characteristics is based on comparison with synthetic signals simulated on models with and without topography.

The comparison between observed and simulated signals has to be carried out very carefully because of the uncertainties of the seismic source and the propagation medium. It is important to emphasize that we do not want to reproduce the recorded signal but rather understand some of its features, such as for example arrival times, waveform complexity, and amplitudes.

For the analysis, a single boulder rockfall is chosen with well separated impacts that can be tracked on video. These criteria are fulfilled by an event that occurred on January 22, 2017, located on the northern crater wall. Figure 13 dust clouds caused by the impacts. The time delay between the actual impact and the visibility of the dust cloud influences the accuracy of the impact time to a similar order of magnitude as the sampling time of 0.5 s between successive snapshots.

The broadband seismic signal of the rockfall shown in Figure 13c is characterized by two main lobes. These lobes are separated by a gap of low seismic energy at around 10:26:32. During this gap, no impact is observed on the video. Thus, the boulder is in free fall before hitting the ground at impact location N2. Afterwards, the rockfall splits into several blocks that continue to move downwards on the debris cone of former rockfalls. At these later times it is very difficult to identify single impacts.

To better distinguish single impacts, the seismic signal is filtered in different frequency bands. Figure 13d The signal filtered in the low frequency band (2-10 Hz) exhibits a smooth amplitude envelope. The two main lobes discussed above can be observed whereas no single pulses can be identified. This signal contains the strongest amplitudes and thus dominates the broadband signal. Short signal pulses emerge in the high-frequency bands. It is clear that seismic sources were already active before impact N1. Unfortunately, they could not be detected on the video. Impacts are possibly hidden behind the clouds on the top of the crater wall. A clear seismic pulse in the frequency range 10-20 Hz can be ascribed to impact N1. It arrives at the station around 0.5 s after the time determined from the video. A second pulse around 1 s later can be ascribed to impact n1. It contains slightly smaller amplitudes. The highest frequency band does not show signals clearly corresponding to these two impacts. This is different for impact N2. Both high-frequency bands show abrupt signal onsets around 1 s after the detection of impact N2 on the video.

The following signal cannot be described as a single pulse but contains several peaks. This raises the question as to whether the source is made up of several impacts or if these peaks result from seismic wave propagation.

Another interesting observation concerns the impact-generated frequencies. As we can see, impact N1 is barely detectable in the highest frequency range (20-40 Hz), whereas impact N2 produces clear signals in both high-frequency bands (10-20 Hz and 20-40 Hz).

Considering the changing source-receiver distance, we would expect the contrary as N1 is slightly closer to station BON than N2. If we assume that the properties of the boulder and of the underlying ground are identical for both impacts, the change in frequency content must be related to the impact velocity. As the boulder accelerates between impact N1 and impact N2, the higher velocity at impact N2 results in a shorter collision time and therefore generates higher frequencies according to Hertz contact theory, introduced below.

Hertz contact theory

To predict relative amplitudes of signals generated by impacts N1 and N2, the respective impact forces of the boulder on the ground are estimated. [START_REF] Farin | Characterization of rockfalls from seismic signal: Insights from laboratory experiments[END_REF] used the theory of [START_REF] Hertz | Die lieber die Berührung fester elastischer Körper[END_REF] to describe the force of an elastic sphere impacting a solid elastic surface. After successfully applying the theory on seismic signals generated in laboratory experiments, they analyzed real-size rockfall experiments carried out by [START_REF] Dewez | OFAI: 3D block tracking experiment on a weathered volcanic rock slope of Tahiti, French Polynesia[END_REF]. Here we estimate the impact forces in a similar fashion, assuming a spherical boulder of radius R and mass m, where m = ρ 4/3 πR 3 with rock density ρ. The maximum impact force F 0 exerted by the sphere perpendicularly to the plane can then be expressed as [START_REF] Johnson | Contact Mechanics[END_REF])

F 0 = 4 3 ER 1/2 δ 3/2 max , (3)
where δ max is the maximum indentation depth manuscript submitted to JGR: Solid Earth with impact speed v n normal to the plane. E is the effective Young's modulus 1/E * = (1 -ν 2 s )/E s + (1 -ν 2 p )/E p , where ν s , ν p , E s , and E p are Poisson's ratio and Young's modulus of a sphere and an impacted plane, respectively.

δ max = 15mv 2 n 16ER 1/2 2/5 , (4)
Concerning the frequency content of the impacts, we analyze the contact duration of the impacts. As proposed by [START_REF] Johnson | Contact Mechanics[END_REF], the temporal evolution of the Hertzian impact force F H can be approximated by

F H (t) ≈ F 0 sin(πt/T c ) 3/2 , 0 ≤ t ≤ T c . (5)
The force-time function and its frequency spectrum are shown in Figure 14 The spectral amplitude decays in the form of a power law above corner frequency [START_REF] Johnson | Contact Mechanics[END_REF] showed that the impact duration can be approximated by means of maximum indentation depth δ max and impact normal speed v n as

f c = 1/T c .
T c ≈ 2.94 δ max v n . (6)
The values of impact parameters for N1 and N2 are estimated assuming that boulder and ground properties are identical between the two impacts. After careful consideration, the uncertainties for all parameters is uniformly set to ±50%, which is found to be reasonable as a maximum threshold from a physical perspective and also allows numerical comparison of the respective contribution of each parameter on the total uncertainty of impact force and corner frequency.

To estimate impact speed v n normal to the slope, a sub-vertical fall of the boulder before collision is assumed for both N1 and N2 with vertical speed v c at the time of collision. Slope angles at the impact positions are inferred from the DEM to be around α = 15 • . The normal impact speed can then be calculated as v n = v c sin α. To estimate v c for N1 and N2, height differences between the impacts are determined from the DEM using the impact positions estimated from the video. As labeled in Figure 13a, we , where g = 10 m s -2 is acceleration due to gravity. Hence, a vertical speed v c,2 = 65 m s -1 is found for impact N2. This leads to normal impact speeds of v n,1 = 4 m s -1 and v n,1 = 17 m s -1 for impact N1 and N2, respectively.

The boulder size is approximated from camera snapshots. The dust cloud caused by impact N2 in Figure 13a has an estimated length of 5 m. As only the dust clouds and not the boulder itself can be seen on the video, the boulder size is assumed to be less than 1 m. As a lower bound, a size bigger than 0.2 m is assumed necessary to generate a seismic signal clearly above the ambient noise level. We therefore estimate the boulder size to be the mid-point between these two limits, i.e. 0.6 m. As a perfect sphere is considered in the calculations, the effective radius is hence estimated to be R = 0.3 m.

The rock density is estimated to be ρ = 2000 kg m -3 , as in the simulations, which results in a boulder mass of m = 226 kg. A typical effective Young's modulus of E = 10 MPa is applied as proposed by [START_REF] Farin | Characterization of rockfalls from seismic signal: Insights from laboratory experiments[END_REF].

The maximum impact force F 0 can now be calculated using equation 3. We find 84 kN and 487 kN for impacts N1 and N2, respectively. The maximum deviations, summarized in Table 2 and ranging from -95% to +450%, are estimated numerically by varying each parameter by ±50%. Figure 15a breaks down the contribution of each parameter to the maximum uncertainty. Note that the relative errors are the same for both impacts N1 and N2. We can observe that a variation of impact speed v n has the greatest effect on the impact force, followed by the rock density ρ. b Parameters for impacts N1 and N2: vertical impact speed v c , angle α between the slope and the vertical, impact speed v n normal to the slope, maximum indentation depth δ max , impact force F 0 (with maximum deviations), contact time T c , and corner frequency f c (with maximum deviations). The impact duration is estimated to be 0.038 s and 0.029 s for N1 and N2, respectively. As f c = 1/T c , it follows that the high-frequency content of the impacts are lim-ited by corner frequency 26 Hz and 35 Hz, respectively, with maximum deviations ranging from -62% to +236%, see Table 2. The contribution from each parameter is indicated in Figure 15b for impacts N1 and N2. In contrast to the impact force, the frequency content is least sensitive to the normal impact speed v n . It is most sensitive to rock density ρ and Young's modulus E.

An important result is that Hertz contact theory predicts a higher frequency content for N2, related to the higher impact velocity. This agrees with the observed waveforms in Figure 13d: impact N1 can hardly be detected in the high-frequency band (20-40 Hz), whereas impact N2 shows a clear pulse despite the slightly bigger source-receiver distance. The theoretical values agree well with the observations, predicting frequencies up to 26 Hz and 35 Hz for N1 and N2, respectively.

Comparison of observed and synthetic waveforms

Previously, we tried to associate pulses in the observed seismic signal to impacts observed in the video by crudely interpreting the signal after the time of impact. We will now use numerical simulations to obtain insights into travel times and expected waveforms. As already mentioned, the intention is not to reproduce observed waveforms, but rather to understand which signal characteristics can be ascribed to a single impact. Analyzing the synthetic seismograms, it can generally be observed that N1 produces smaller amplitudes compared to N2. This is due to the estimated impact forces of 84 kN and 485 kN, respectively (see Table 2).

While seismograms from the model with a flat surface keep approximately the same relative amplitudes between N1 and N2 at all different stations, seismograms from the model with topography show more variability. For example, at station BON, the amplitudes of impact N1 are much bigger for the model with topography than for the flat model. This corresponds better to the real observations, where the maximum amplitude of impact N1 is comparable to the maximum amplitude of impact N2. In contrast, at station DSO, the signal of impact N1 is very small as opposed to the signal of N2 for the model with topography. Again, this corresponds well with the observations. As impacts N1 and N2 are located very close to each other, the relative amplitudes on the flat model are mainly determined by the relative impact force. In contrast, surface topography influences both the relative (vertical) source position and the propagation path, so that small source displacements can cause local amplification or deamplification of the signal. Measuring the signal at station DSO, the source moves from a deamplified zone at N1 towards an amplified zone at N2, as suggested by the reciprocal amplification pattern provided in the Supporting Information.

From the simulations with a flat surface, three wave packets can be observed following each impact, which are well separated from each other on the more distant stations BOR, DSO, and SNE. These three wave packets correspond to a body wave, a first- A1.

mode Rayleigh wave, and a fundamental mode Rayleigh wave (see e.g. Figure 3). The arrival time of the first-mode Rayleigh wave is in good agreement with the first major pulse after each impact. This suggests that the Lesage generic velocity model represents the shallow subsurface velocity around Dolomieu crater reasonably well. However, for the flat model, the amplitude of the first-mode Rayleigh wave is consistently smaller than the amplitude of the fundamental mode. A corresponding amplitude variation cannot be identified on the real signals. In contrast, simulations on the model with topography generate more complex waveforms. This increased complexity corresponds better to the observed signals, even if the waveforms do not fit perfectly. The variation of the force direction modifies the waveforms more than in the flat case. Also, waveforms vary greatly manuscript submitted to JGR: Solid Earth from station to station. This is not observed in the flat case, in which the waveforms are very similar for stations at comparable source-receiver distances (i.e. BOR, DSO, and SNE). Similar observations are made for synthetic signals generated by neighboring source positions (10 m spacing), exhibiting much more variability for the case of topography, see Supporting Information.

It cannot be excluded that the observed amplitude variations and waveform complexity can be explained by heterogeneous velocity structures or the superposition of multiple (here undetected) sources. Nonetheless, the simulations indisputably show that these effects are caused by surface topography, for which high-resolution models are available.

Consequently, surface topography must be taken into account when interpreting highfrequency rockfall seismic signals and multiple signal pulses must not directly be interpreted to be caused by multiple impacts.

Finally, observed and synthetic seismograms are compared without normalization.

In this way, the absolute signal amplitudes calibrated by the Hertz impact force are evaluated. The total value of the acting force as well as its direction are determined by a vector sum of the forces normal and tangential to the slope. Tangential force F t is inferred from the maximum normal impact force

F n = F 0 assuming Coulomb friction F t = µF 0 ,
where µ is the material-specific friction coefficient. We assign µ = 0.7, a typical value used for rockfalls at Dolomieu crater (e.g. [START_REF] Hibert | Automated identification, location, and volume estimation of rockfalls at Piton de la Fournaise volcano[END_REF]. The resulting signal amplitudes for model simulations with a flat surface and with topography are compared with the observed rockfall signals in Figure 17. A1.

F n + F t F n + F t flat topo BON F n + F t F n + F t flat
Despite the uncertainties on impact force F 0 and the assumptions on elastic contact and Coulomb friction, which are certainly violated by a certain amount of plasticity and surface roughness, we are able to achieve order-of-magnitude agreements with the observed signals. The application of the Hertz contact theory to calibrate the simulations therefore offers a solution to this poorly constrained problem.

manuscript submitted to JGR: Solid Earth

When comparing the amplitudes between the model with a flat surface and the model with topography, one notices that the relative amplitude changes between the stations.

While amplitude decreases linearly with the distance from the source for the flat model, this is not the case for simulations with topography. For example, for the flat model, the maximum amplitude decreases drastically from station BON to station DSO. In contrast, for the model with topography, the signal amplitude of impact N2 at station DSO is only about half that of station BON, which can be assigned to topographic amplification.

Conclusion

Using Spectral-Element-based simulations with high-resolution surface topography, we were able to broadly explore the influence of topography on the propagation of wave fields generated by surface point loads and to achieve order-of-magnitude agreement with seismic measurements from rockfalls at Dolomieu crater.

In the first half of the present study we numerically simulated the topography effect for different earth models. It was shown that the topography effect is significantly altered depending on the underlying velocity model. For example, for a homogeneous model, seismic energy is easily directed downwards by the topography into the subsurface. This is contrary to a velocity-depth profile with a strong gradient, as proposed by [START_REF] Lesage | A generic model for the shallow velocity structure of volcanoes[END_REF] for the shallow velocity structure of volcanoes, for which more seismic energy is trapped close to the surface leading to more complex spatial distributions of amplification.

Studying the geometric features of Dolomieu-like crater topographies suggests that curvature variations affect the seismic wave field more than depth variations. However, the locally inherent complexity of the topography prohibits a holistic generalization of its effect on the seismic wave field. Nevertheless, our approach provides a methodology to quantify the influences of the topography. This methodology can be used to other study sites by changing the model parameters regarding the domain surface and the medium properties, which are inevitably coupled. In this respect, scattering from 3D soil heterogeneities, which can strongly disturb wave propagation, must be considered. Regarding the study site at Dolomieu crater, previous studies and our findings indicate that scattering does not dominate the signals.

Local site amplification at the stations could be taken into account by using amplification factors estimated from VT events, which significantly improved the agreement between simulations and observations. In general, however, it cannot be guaranteed that surface waves will receive the same amplification as a vertically incident wave field and it is recommended that wave propagation, including subsurface heterogeneities, be modeled when the necessary information is available.

In the second half of the present study, the simulations were compared with measurements from rockfalls at Dolomieu crater. Different rockfalls with similar impact locations were investigated by calculating spectral ratios between the stations, thereby extracting the signature of the seismic source. The agreement between the observed spectral ratios indicates identical path effects that were reproduced by simulations using the model with topography. It was further shown that the spectral ratios are dominated by propagation along the topography rather than by the direction of the seismic source. As the latter is hard to estimate precisely, this finding can have practical applications, for example for the locating of rockfalls.

Single impacts were shown to be able to generate complex waveforms with multiple pulses depending on medium properties and topography. For topography, variations of source direction and position can strongly modulate the waveforms.

Using Hertz contact theory, signal features were linked to impact parameters. Estimating the maximum impact force helped to calibrate the simulation and achieve orderof-magnitude agreement with observed signal amplitudes. Associating higher frequencies with increased impact speed explained the observed frequency content of rockfall impacts. lected from a catalog compiled by [START_REF] Duputel | Constraining Spatiotemporal Characteristics of Magma Migration at Piton De La Fournaise Volcano From Pre-eruptive Seismicity[END_REF]. To compute the amplification manuscript submitted to JGR: Solid Earth tral amplitudes of both VT signals and H/V noise ratios. The spectral ratios are computed from FFT spectra after applying the smoothing function proposed by [START_REF] Konno | Ground-motion characteristics estimated from spectral ratio between horizontal and vertical components of microtremor[END_REF] to avoid spurious fluctuations. Figure C1b shows the mean spectral amplification functions and their standard deviation calculated from all VT events for all components.

Strongest amplification is experienced by station SNE with factors up to 7 on its horizontal components. This can explain the strong mismatch between observations and simulations that are visible in Figure C1a. The vertical component of single-component station DSO also seems to be amplified with a peak around 9 Hz. Less evidence of amplification is found for station BOR, except for its north-component which is amplified by a factor of 2 for frequencies above 5 Hz.

Appendix D Observed and simulated spectral ratios for rockfall sources in the southwest

To reinforce the findings of section 5.2.2 that spectral ratios are characteristic of the source position and can be reproduced when surface topography is taken into account, the same analysis is carried out here for rockfall sources located on the southwestern crater wall. Snapshots taken from camera CBOC of the three observed rockfalls are shown in The camera images reveal that each of the rockfalls involves at least two boulders moving downslope simultaneously. While the boulders of rockfall 1 originate from below the camera position, boulders of rockfall 2 and 3 come from the right-hand border of the image. At the time of the snapshot, the trajectories of the three rockfalls cross.

From a window of ±4 s around this time, the spectral ratios are computed from the observed signals and shown in Figure D2 (blue, site-effect corrected). As for the rockfalls analyzed above, the spectral ratios from the three events in the southwest are very similar to each other across the whole frequency range for all station pairs. The spectral ratios are now compared to simulations using the model with Dolomieu topography. As above, three input force configurations are tested (i.e. a vertical force, a force normal to the slope and a force tangential to the slope in direction of the strongest gradient) to investigate the dependency of the ratios on the source direction.

We can generally observe that the simulated spectral ratios agree very well with the observed spectral ratios. Changing the source direction does not essentially influence the spectral ratios, except for frequencies below 5 Hz, which is similar to the observations in Figure 12. The similarity at higher frequencies suggests that the ratios are dominated by propagation along the topography rather than by the source mechanism. The strongest deviation between observations and simulations at high frequencies is visible on ratio SNE/BON for the east-component. In comparison with the observed spectral ratios, the simulated amplitudes measured at station SNE are strongly underestimated with respect to station BON. This could be caused either by soil heterogeneities, which are not considered in the simulations, or by uncertainties of the source position and the fact that the rockfall contains at least two boulders which simultaneously impact the ground.

Analysis of the rockfalls located in the southwestern part of Dolomieu crater supports the findings of section 5.2.2 indicating that the spectral ratios are characteristic of the source location, and can be reproduced by taking into account the surface topography, while source direction is not dominant, in particular at high frequencies.

  ://doi.org/10.1002/essoar.10502632.2 | CC_BY_4.0 | First posted online: Wed, 14 Oct 2020 00:47:55 | This content has not been peer reviewed.

FigureFigure 1

 1 Figure2ashows a cross-section through the spectral-element mesh with Dolomieu surface topography. The dimensions of the domain measure x = 2100 m (easting), y =

Figure 2

 2 Figure 2. a) Cross-section of the SEM mesh through Dolomieu crater with a topography resolution of 10 m. Perspective as seen from the East with Bory crater located in the background. The color map corresponds to the Lesage generic velocity model (see section 3.2). The buffer layer 100 m below the surface dampens small-scale topography variations. The zone of refinement at 150 m below the surface connects elements with 10 m and 30 m side lengths. 160 m wide PML boundaries are attached to the sides and bottom of the domain. b) S-and P-wave velocity depth profiles for the (1) homogeneous model (v S,1 and v P,1), (2) model with shallow S-wave velocity layer (v S,2 and v P,2), and (3) Lesage generic velocity model (v S,3 and v P,3). The shaded zone (v S,M o) is extracted from the inverted 3D S-wave model of Mordret et al. (2015). c) Theoretical dispersion curves of the Lesage generic model for the fundamental (R0) and first-mode (R1) Rayleigh wave velocity together with picked dispersion curves from a mini-array around station BON. The errors are estimated from the uncertainty during dispersion curve picking.

Figure 3 .

 3 Figure 3. Influence of topography resolution on synthetic seismograms from the Lesage generic velocity model, recorded at stations BON, BOR, DSO and SNE. a) Comparison of syn-

Figure 4

 4 Figure 4 shows synthetic seismograms recorded on the surface along an array crossing the source position, Dolomieu crater and station BON (see inset for location of the array). Snapshots of the propagating seismic wave field on a cross-section along the array are shown below the seismograms. All amplitudes correspond to vertical ground velocity. In order to enhance visibility of the wave field over time, the simulations here are carried out without intrinsic attenuation.

Figure 4 .

 4 Figure 4. Wave propagation from a vertical surface load for different velocity models. Synthetic seismograms (top row ) recorded at an array crossing the source, Dolomieu crater and station BON (see inset) for the (1) homogeneous model (left), (2) model with a shallow S-wave velocity layer (middle), and (3) Lesage generic velocity model (right). The seismic traces are normalized with respect to themselves and show vertical ground velocity. Snapshots of the wave field on cross-sections along the same array are shown below the seismograms, corresponding to the times marked by red dashed lines. Annotations denote P-wave (P), S-wave (S), P to S converted wave (PS), Rayleigh to S converted wave (RS), Rayleigh wave (R), reflected Rayleigh wave (Rr), diffracted Rayleigh wave (Rd), and diagonally traveled Rayleigh wave (Rf). SEM configurations correspond to 3, 6 and 11 in TableA1. The absence of intrinsic attenuation for enhanced

Figure 5 .

 5 Figure 5. Topographic amplification from for a vertical point force. Amplification for of vertical PGV (top) and total kinetic energy (bottom) is calculated relative to a flat reference model for the homogeneous model (left), the model with a shallow low velocitylow-velocity layer (middle) and the Lesage generic velocity model (right). SEM configurations correspond to 1, 2; 4, 5;

  Figure 6 illustrates vertical PGV amplification (left) and energy amplification (right) for a wave field generated by a horizontal surface force polarized in the north direction.A strong directionality is visible in the PGV amplification pattern. This is because for the flat reference model, a horizontal source does not generate vertical seismic energy perpendicular to its direction. Topography however can change this by conversion from transverse energy or diffraction of wave paths.

FFigure 6 .

 6 Figure 6. Topographic amplification for a horizontal point force in the north-direction. Amplification of vertical PGV (left) and total kinetic energy (right) is calculated relative to a flat reference model for the model with the Lesage generic velocity profile. SEM configurations correspond to 9 and 10 in Table A1. The black arrow illustrates the source position and its direction. Green triangles mark station locations. Annotations give ratios measured at the station locations as well as the percentage of topographic amplification. Neighboring contour lines differ by 60 m in elevation.

Figure

  Figure 7a compares synthetic seismograms recorded along arrays on the domains with a flat surface, a planar rough surface, and crater topography.For the model with the flat domain, we can identify dispersive fundamental and first-mode Rayleigh waves as well as body waves. Introducing surface roughness leads to strong scattering and hence prolonged ground shaking. The two Rayleigh modes are no longer clearly separated, even though the propagation of the main energy from the fundamental mode can be identified. Introducing the crater topography adds more complexity. In particular the steep crater walls distort the propagating wave field, as already observed for the real crater topography (see Fig.4).

Figure 7

 7 Figure 7. a) Synthetic seismograms of vertical ground velocity from models with a flat surface (left), a rough planar surface (middle), and a synthetic crater shape (right). Seismograms are normalized and recorded along the surface profiles illustrated by the blue curves on top. The yellow star marks the position of the vertical source. The spurious reverberations observed for the flat surface model after the signal (> 6 s) are trimmed for the analyses. b) Energy amplification in different frequency bands for the model with a rough surface (left) and with a synthetic crater (right) relative to the flat reference model for frequency bands 3-7 Hz and 13-17 Hz. The arc-like features at 3-7 Hz in the top corners are numerical artefacts. SEM configurations correspond to 12, 13 and 14 in TableA1.

1 )Figure 8 .

 18 Figure 8. Left: Profiles through the synthetic crater topographies. a) Crater depths vary by ±0.3 from small to big with fixed curvature. b) Curvatures vary by ±0.3 from weak to strong with fixed crater depth. Red dashed lines correspond to a profile through Dolomieu crater and its corresponding curvature. Right: Comparison between energy amplification for crater geometries with smallest and biggest depths as well as weakest and strongest curvature. Contour lines mark elevation differences of 50 m and the yellow star denotes the source. Note that spurious blue dots inside the crater (especially at steep flanks for the big depth) were caused by numerical measurement problems at these positions. SEM configurations correspond to 12 and 14 in Table

  ://doi.org/10.1002/essoar.10502632.2 | CC_BY_4.0 | First posted online: Wed, 14 Oct 2020 00:47:55 | This content has not been peer reviewed.

Figure 9 .

 9 Figure 9. Three similar rockfalls on the southern wall of Dolomieu crater, corresponding to rockfall location 1 in Figure 1b. The events occurred on: 1) February 28, 2016 at around 11:47, 2) February 28, 2016 at around 12:46, and 3) February 18, 2016 at around 12:27. a) The total trajectory of each event (seen from camera SFRC). The approximate starting positions at the top of the crater wall are indicated by white arrows. b) Snapshots (seen from camera SFRC) at a chosen time at which all three rockfalls are at comparable positions. c) The rockfall seismic signals, the red dotted lines indicating the times of the snapshots. Time windows R1, R2, and R3(blue-shaded zones) are defined ±4 s around these times. The corresponding locations of these time windows are also indicated as blue-shaded zones on the trajectories. The same holds for reference time window C1 (magenta-shaded zone), which corresponds to the beginning of event 1.Noise time window N is taken from recordings before event 1.

Figure 10 .

 10 Figure 10. Spectral ratios from rockfall seismic signals recorded at stations BOR (3 components), DSO (1 component) and SNE (3 components) relative to station BON for vertical-(top), north-(middle) and east-(bottom) components. Time windows TW-R1, -R2, and -R3 correspond to rockfalls 1, 2, and 3 as defined in Figure 9. Time windows TW-N and TW-C1 correspond to noise recordings and the beginning of rockfall 1, respectively. Time window TW-C2 is taken from a rockfall on the southwestern crater wall, corresponding to rockfall location 2 in Figure 1b.

  ://doi.org/10.1002/essoar.10502632.2 | CC_BY_4.0 | First posted online: Wed, 14 Oct 2020 00:47:55 | This content has not been peer reviewed. manuscript submitted to JGR: Solid Earth ulations. Therefore, in order to enhance comparability between the observed and simulated spectral ratios, the recorded signals are corrected using site amplification factors estimated from volcano-tectonic (VT) signals. The spectral amplification curves are calculated and discussed in Appendix C where we also show a comparison between simulated and uncorrected observed spectral ratios. The observed inter-station ratios presented here are corrected by deconvolution of the recorded signals with the corresponding amplification factors.

Figure

  Figure 11c compares spectral ratios BOR/BON from the observed signals with simulated ratios of differently polarized sources on a model with a flat surface. The source directions are determined from the slope of Dolomieu topography at the corresponding position before implementation on the flat domain. For spectral ratios of the vertical component (left), a tangential force direction results in much smaller values compared to the other sources. As the slope dips northwards, the tangential force is orientated in the north-

Figure 11

 11 Figure 11. a) Seismic radiation patterns from a vertical source (top) and from a horizontal force (bottom) for ground velocity vZ , vN , and vE of vertical-, north-, and east-component, respectively (red for positive and blue for negative amplitudes). b) Forces generated by a rockfall impact. The red dotted line illustrates the trajectory of a bouncing boulder. The impact generates force Fn normal to the slope. Depending on the boulder velocity tangential to the slope and on the friction coefficient µ, a tangential force Ft = µFn is generated (assuming Coulomb friction). Normal and tangential forces add up to resulting force Fr. c) Comparison of spectral ratios BOR/BON from real signals TW-R1, -R2, -R3 (as in Fig. 10, site-effect corrected) and from simulations on the flat domain with varying source direction: vertical force, normal force and tangential force according to the Dolomieu topography at the corresponding position. The shaded zones of the simulated ratios indicate the standard deviation around the mean value from seven neighbouring source positions (Fig. 1b, picked sources). SEM configurations correspond to 15 in TableA1.

Figure 12 .

 12 Figure 12. Spectral ratios BOR/BON, DSO/BON and SNE/BON calculated from real signals TW-R1, -R2, -R3 (as in Fig. 10, site-effect corrected) and from simulations on the domain with Dolomieu topography for the vertical (top), the north (middle), and the east (bottom) component. Simulations are carried out with varying source directions: vertical force, force normal to the slope and force tangential to the slope. The shaded zones of the simulated ratios indicate the standard deviation around the mean value from seven neighbouring source positions (Fig. 1b, picked sources). SEM configurations correspond to 16 in TableA1.

Figure 13 .

 13 Figure 13. Single boulder rockfall on January 22, 2017. a) Camera snapshot taken shortly after impact N2. Estimated vertical distance between impacts and estimated slope angle to the vertical at impact positions. b) Location of impacts N1 and N2 in Dolomieu crater. c) Vertical ground velocity recorded at closest station BON in frequency band 2-40 Hz. The red-shaded area illustrates the time window of the graph below. Dashed lines mark impact times N1 and N2

  compares signals band-pass filtered at 2-10 Hz, 10-20 Hz, and 20-40 Hz. The relative scales of the normalized signals can be inferred from the gray bars plotted at the beginning of the signal as well as from the spectrogram below.

  ://doi.org/10.1002/essoar.10502632.2 | CC_BY_4.0 | First posted online: Wed, 14 Oct 2020 00:47:55 | This content has not been peer reviewed.

Figure 14 .

 14 Figure 14. Hertzian impact force and corresponding frequency spectrum. Left: Hertzian force-time function FH normalized by maximum impact force F0 versus impact duration Tc, which represents the time during which the two bodies are in contact. Right: Frequency spectrum of the force-time function. The inverse impact time 1/Tc is related to the corner frequency fc after which the spectral amplitude decays in the form of a power law.

  find a height difference of around H 1 = 15 m between N1 and n1 and a height difference of around H 2 = 140 m between N1 and N2. Impacts N1 and n1 are detected 1 s apart. Assuming an approximately constant velocity during this short time window, the vertical speed for impact N1 is v c,1 = 15 m s -1 . For impact N2, acceleration during the long free fall cannot be neglected. The speed is thus derived by v c,2 = v c,1 +(2g(H 2 -manuscript submitted to JGR: Solid Earth H 1 ))0.5

Figure 15 .

 15 Figure 15. Influence of individual parameters on uncertainty of impact force F0 (a) and corner frequency fc (b). Impact parameters include boulder radius R, rock density ρ, normal impact speed vn and Young's modulus E. Uncertainties are estimated numerically by varying each parameter by ±50%. Note that the relative errors for N1 and N2 are identical.

  Observed and synthetic signals are compared in the frequency band of 10-20 Hz, in which we identify short signal pulses caused by the rockfall impacts. At the same time, 20 Hz constitutes the upper frequency limit of our simulations. In order to ensure comparability, the observed signals are corrected with the site amplification functions calculated in Appendix C and subsequently convolved with the 7 Hz Ricker wavelet used in the simulations. The simulated signals are calibrated with the maximum impact forces for N1 and N2 estimated above using Hertz theory.Initially, observed and synthetic seismograms are normalized for easier comparison of the waveforms. Figure16compares vertical ground velocity recorded at stations BON, BOR, DSO, and SNE with simulations from impacts N1 and N2 for models with a flat surface and with Dolomieu topography. Source positions of the two impacts are estimated from the videos, see Figure13b. As the exact source direction of the real impacts are unknown, the variability of synthetic waveforms is shown for different force directions, namely a vertical force F z , a force F n normal to the slope and a force F t tangential to the slope.

Figure 16 .

 16 Figure 16. Comparison of recorded rockfall signal (blue, for a time window as in Figure 13d) with synthetic waveforms of impacts N1 and N2 for the model with a flat surface (green) and with Dolomieu topography (orange). The variability of the synthetic waveforms is demonstrated for a vertical force Fz, a force Fn normal to the slope and a force Ft tangential to the slope. Red vertical lines indicate impact times N1 and N2 from the video. Corresponding source positions are shown in Figure 13b. All signals show vertical ground velocity and are normalized by their maximum. SEM configurations correspond to 15 and 16 in TableA1.

FFigure 17 .

 17 Figure 17. Amplitude comparison of recorded rockfall signal (blue, for a time window as in Figure 13d) with synthetic waveforms of impacts N1 and N2 for the model with a flat surface (green) and with Dolomieu topography (orange). The synthetic seismic source is constructed by summing force Fn and Ft normal and tangential to the slope. Fn corresponds to the maximum impact force F0 as shown in Table 2. Red vertical lines indicate impact times N1 and N2 from the video. Corresponding source positions are shown in Figure 13b. All signals show vertical ground velocity and are normalized by their maximum. SEM configurations correspond to 15 and 16 in TableA1.

  Figure C1. a) Comparison of spectral ratios as in Figure 12 but with before site effect correction of the recorded signals (blue). This leads to partly strong deviations between observations and simulations. b) Spectral amplification functions relative to reference BON for ground velocity of vertical (top), north (bottom left) and east (bottom right) component. The blue-shaded zone indicates the standard deviation of the amplification distribution from all VTs.

  957 functions, BON was qualified as an adequate reference station based on the low spec-958 -34-ESSOAr | https://doi.org/10.1002/essoar.10502632.2 | CC_BY_4.0 | First posted online: Wed, 14 Oct 2020 00:47:55 | This content has not been peer reviewed.

Figure

  Figure D1, together with the generated seismic signals recorded on the vertical components. For the times of the shown images, marked on the seismic signals by the vertical dotted red lines R1, R2, and R3, all the rockfalls are located in the same area.

Figure D1 .Figure D2 .

 D1D2 Figure D1. Top: Snapshots taken from camera CBOC of three rockfalls times for which all the rockfalls are in comparable locations. Positions and directions of the boulders are indicated by red circles and arrows. The trajectory of the rockfall on 13 December, 2016 is indicated as event 2 on the map in Figure 1b Bottom: Corresponding seismic signals (vertical velocity). The vertical dotted lines R1, R2 and R3 mark the time of the camera snapshot shown above. The blue-shaded zones display the time windows of ±4 s around R1, R2, and R3 in which spectral station ratios of the signals are computed.

Table 1 .

 1 Parameters of velocity models for the SEM simulations a

	plemented in Geopsy software (www.geopsy.org). Theoretical dispersion curves are cal-
	culated from the Lesage generic model using modal summation from Computer Programs
	in Seismology (Herrmann, 2013). The measured values agree well with the fundamen-
	tal Rayleigh velocity dispersion curve. No coherent dispersion curves could be picked above
	6 Hz because of the minimum mini-array aperture of 30 m.
	Despite missing measurements above 6 Hz, the Lesage generic model is assumed
	to be the most reasonable model for the shallow high-frequency velocity structure of Piton
	de la Fournaise volcano because it is based upon measurements at comparable volcanoes.

Table 2 .

 2 Hertz impact parameters b

	v c	α	v n	δ max	F 0		T c	f c
	N1 15 m s -1 15 •	4 m s -1 0.05 m	84 +376 -79	kN	0.038 s	26 +61 -16 Hz
	N2 65 m s -1 15 • 17 m s -1 0.16 m	487 +2186 -460	kN 0.029 s	35 +82 -21 Hz

  ESSOAr | https://doi.org/10.1002/essoar.10502632.2 | CC_BY_4.0 | First posted online: Wed, 14 Oct 2020 00:47:55 | This content has not been peer reviewed.

		manuscript submitted to JGR: Solid Earth
	N1	N2		N1	N2
	F t			F t
	F n			F n
	F z			F z
	F t			F t
	F n			F n
	F z			F z
	N1	N2
	F t
	F n
	F z
	F t
	F n
	F z
	10:26:28	10:26:32	10:26:36	10:26:28	10:26:32	10:26:36
		Time (s)			Time (s)
			-28-

estimated from the video. d) Comparison of frequency bands 2-10 Hz, 10-20 Hz, and 20-40 Hz. Signals are normalized to their maximum and the gray bars on the left indicate the relative scaling. e) Time-frequency representation of the rockfall signal (calculated using the Stockwell transform).

der impacts, N1 and N2, around 4 s apart, are analyzed. A minor impact n1 is observed 1 s after impact N1. It will be used later to estimate the fall velocity of the boulder. Note that the impact times are estimated from the video according to the appearance of small -24-ESSOAr | https://doi.org/10.1002/essoar.10502632.2 | CC_BY_4.0 | First posted online: Wed, 14 Oct 2020 00:47:55 | This content has not been peer reviewed. manuscript submitted to JGR: Solid Earth

Acknowledgments

We would like to thank Jean-Pierre Vilotte and Hugo Martin for fruitful ideas. We are very grateful to Victor Tsai and an anonymous reviewer for their critical and constructive remarks which substantially contributed to improve the quality of the paper. We thank the whole team at the OVPF that provided the excellent field data for this study. The seismic data used in this study, were acquired by the Volcanological and Seismological Observatory of Piton de la Fournaise (OVPF)/Institut de Physique du Globe de Paris (IPGP) via the VOLOBSIS Portal: http://volobsis.ipgp.fr/query.php. Camera data are stored in http://doi.org/10.5281/zenodo.4031816. Data from the simulations are available from http://doi.org/10.5281/zenodo.3949826. Numerical computations were partly performed at S-CAPAD (Service de calcul parallèle et de traitement de données en sciences de la Terre), IPGP, France, as well as at CCIPL (Centre de Calcul Intensif des Pays de la Loire), Université de Nantes, France. This work was funded by the ERC Contract No. ERC-CG-2013-PE10-617472 SLIDEQUAKES.

manuscript submitted to JGR: Solid Earth

The combination of Hertz contact theory and wave propagation simulations is an important step for the interpretation of rockfall seismic signals based on the underlying physical processes. The Hertz impact theory is frequently used to predict the impact force of rockfalls, for example for the design of protective structures (e.g. [START_REF] Volkwein | Rockfall characterisation and structural protection -A review[END_REF]. Also, laboratory experiments show the validity of Hertz theory concerning the waves generated by the collision of a ball on a massive plate (e.g. [START_REF] Mclaskey | Hertzian impact: Experimental study of the force pulse and resulting stress waves[END_REF] or grains on a plate [START_REF] Bachelet | Elastic wave generated by granular impact on rough and erodible surfaces[END_REF]. However, only a few studies apply the theory to seismic signals from real-scale rockfalls [START_REF] Farin | Characterization of rockfalls from seismic signal: Insights from laboratory experiments[END_REF][START_REF] Hibert | Single-block rockfall dynamics inferred from seismic signal analysis[END_REF].

A limiting factor is the complexity of the rockfall source that often consists of multiple simultaneous impacts. For this reason, application to artificially triggered rockfalls, ensuring separate impacts of a single boulder ,would help validate the Hertz theory in the field and enhance our understanding of real impact processes. events which are centered around 2 km below Dolomieu crater. Thirty-six events are se-

Appendix A Configuration of SEM simulations

Appendix B Energy amplification in different frequency bands