Unraveling Information from Seismic Signals Generated by Gravitational Mass Movements

Abstract

A practical analysis of the spectrograms of the seismic data generated by gravitational mass movements (GMMs), such as snow avalanches, landslides, lahars, and debris flows recorded on one sensor, is presented. The seismic signal produced by these movements is analyzed in terms of the shape of the initial section of the spectrogram, which corresponds to the start of the movement of the gravitational mass. The shape of the envelope of the spectrogram is a consequence of the progressive reception of high-frequency energy in the signal as the gravitational mass (GM) approaches the sensor because of the attenuation properties of the seismic waves in the ground. An exponential law was used to fit this envelope of the onset signal. The proposed methodology allows us to obtain the propagation characteristics of different types of GMM. The analysis of the adjusted parameters for different types of GMM allows us to assert that differences of one order of magnitude exist in the values of these parameters depending on the type of event. In addition, differences in the values of the exponent were obtained between the events of each type of the analyzed GMM. We present a template of different curves for each type of GMM with the corresponding parameter values that can help professionals characterize a GMM with only one seismic record (one seismic sensor) whenever the mass movement approaches the recording sensor or passes over it.

1. Introduction

Earth surface gravitational mass movements (GMMs), such as landslides, snow avalanches, debris flows, lahars, pyroclastic flows, and rockfalls present economic and social hazards. Understanding and controlling the behaviors of GMMs is crucial to reduce their associated risks. The monitoring of these events helps in the understanding of hazardous areas basically in two ways: (a) to identify the number of events and areas of occurrence, the inputs necessary for hazard assessment, land use planning, etc., and (b) to obtain a better understanding of the characteristics of the mass movement descending a slope to obtain boundary conditions for the control of the models and thus gain further insight into its evolution. A challenge for professionals involved in risk management is the need to perform remote detection and characterization of the GMMs in which they are involved. The use of seismic signals generated by GMMs is an emerging field to control these phenomena. This is a noninvasive method that supplies information on the evolution of phenomena and detects them.

The wave field generated by a GMM recorded at a seismic station was first investigated by [1,2] for surface waves and body waves excited by landslides generated by the 1980 Mount Saint Hellen eruption. Since then, various theoretical and experimental studies involving different types of GMM (e.g., [3,4,5,6,7]) have been carried out. Most of these studies use more than one seismic sensor and consider and analyze only the generation of the wave field at the source. However, they do not consider, when observing the wave field at a distance, the effect of the ground on the wave propagation and the effect of the site (or site response) on sensor location. These two effects are very well known in the field of seismic engineering (e.g., [8] and references therein). The propagation effect is related to wave attenuation by geometrical spreading and intrinsic attenuation. The consequence of the first effect is the decrease in the signal amplitudes with distance. The second effect produces, with increasing distance, a decrease of the high-frequency content of the signal. As a result, the recorded seismic waves may have different amplitudes and frequency contents than those originating from the source. The consequences of the site effect are not as direct and affect the wave field at the recording site. The topography, ground heterogeneities, and distribution of the ground layers near the sensor can cause amplification, or the opposite, of the amplitudes of different frequencies because of wave interference (constructive or destructive). These effects, if not considered, can cause problems when the characteristics of GMMs are deduced because they can disguise the characteristics at the source of the studied phenomenon.

The evolution of the amplitude of the time series (e.g., [6,9]), both of high-frequency and low-frequency content, is commonly used to obtain information on the dynamics of GMMs (rock falls, landslides, lahars, rock avalanches, and submarine slumps). In [10], a revision of the state-of-the-art of the use of seismic signals related to mass movements associated with volcanoes, which is also valid for other GMMs, can be found. The application of seismic methods continues, and [11,12] are examples of debris flow and flood detection, respectively. Recently, distributed acoustic sensing (DAS), referred to as photonic seismology [13], has been used for the real-time seismic location of GMMs (e.g., [14] for rockfalls or [15] for snow avalanches); however, there is still much to be resolved to equate it with the advantages of classical sensor distributions for exploring GMMs (amplitude control, tractability, mobility, availability, etc., e.g., [15,16]). This is essentially because of the difference in the types of measurement (ground strain rate (DAS) vs. ground particle velocity (seismic)).

In addition to the use of time series (seismograms), the power spectral density (PSD) or spectrum of the signal is also used, but information regarding the frequency evolution over time is lost. Considering the PSD of specific time windows has been proposed for isolating different characteristics of a debris flow [11]. However, little attention has been given to the continuous evolution of the frequency content through the spectrogram. This is a visual representation of the frequency content of a signal as it changes over time; it is also known as a short-time Fourier transform (STFT) and shows the matrix of spectral amplitude values (PSDs) of different correlative time windows. We use the word “spectrogram” instead of “spectrogram of the seismic signals generated by a mass movement” hereafter.

The evolution of the frequency content of a signal over time is a characteristic of each phenomenon. In the case of an earthquake or the impact of a rock, all the frequencies appear simultaneously at the beginning of the spectrogram [17,18], unlike a GMM. An example of the different appearance of the seismogram and spectrogram of an earthquake and a snow avalanche is shown in Figure S1. In these cases, the shape of the initial section of the spectrograms is characterized by an increase in the amplitudes of high-frequency content over time. This shows an envelope that increases exponentially as the avalanche approaches the seismic sensor. In the literature, different examples exist. The spectrogram of the seismic signal produced by a vehicle shows this shape in ([19], Figure 1). This phenomenon is also observed for pyroclastic flows at the Soufrière Hills Volcano in Montserrat in ([4], Figure 11) of and ([20], Figure 6). For snow avalanches, ([21], Figure 8) is a good example. However, in these contributions, the increase in the amplitudes of high-frequency content over time is not mentioned. Nevertheless, in [22], different time window spectra of seismic signals of snow avalanches were used and indicated that the dominant frequencies increase with the approach of the avalanche (to the sensor). This increase in amplitudes in avalanche spectrograms, relating them to the forward motion of the avalanche, was first highlighted in [23].

This paper has two objectives: (a) to demonstrate that all GMMs studied (landslides, avalanches, debris flows, and lahars) exhibit behavior similar for the exponential shape of the envelope of the spectrograms, where their parameters depend on the type of phenomenon and (b) to demonstrate that, as occurs in the case of snow avalanches [24], the variations in these parameters within the same type of GMMs are associated with their evolution. This analysis allowed the development of a template.

This paper is divided into six sections. First, the-state-of-the-art of the use of the spectrogram of seismic signals generated by the GMM is briefly presented. Section 2 introduces the concept of the SON_spectrogram section, the exponential law that governs it, and our objectives. Section 3 provides a description of the data used, a brief explanation of the experimental sites and the characteristics of the different GMMs used. The theoretical approach, a physical approach of the SON_spectrogram interpretation related to the attenuation properties of wave propagation in the ground, is developed in Section 4. The results of the application of the method to different signals generated by snow avalanches at two different sites, lahars, debris flows, and one landslide are presented in Section 5. A template created to classify events and an application to different examples are also presented. The conclusions are provided in Section 6. The acronyms used are indicated at the end of the paper.

2. Outline

In [24], a method to estimate the parameters that define the exponential shape of the SON_spectrogram section was presented. This shape follows an exponential law ($F\left(t\right)=\mathsf{{\rm K}}{e}^{\beta t}),\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{v}\mathrm{a}\mathrm{l}\mathrm{u}\mathrm{e}\mathrm{s}\mathrm{o}\mathrm{f}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{p}\mathrm{a}\mathrm{r}\mathrm{a}\mathrm{m}\mathrm{e}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{s}\mathsf{{\rm K}}\mathrm{a}\mathrm{n}\mathrm{d}\mathsf{\beta }$ change according to the type of snow avalanche. Furthermore, these parameters allowed us to estimate the average speeds of the avalanches that approached the recording locations. In Figure 1, an example of the seismic signal and spectrogram produced by a snow avalanche is presented. Two sections can be distinguished: The initial section of the spectrograms (hereafter referred to as the signal onset (SON_spectrogram)) is characterized by an increase in the amplitudes of high-frequency content over time. After it follows the signal body (SBO) section, where the avalanche passes over or near the sensor, and all the frequencies present maximum amplitudes [24,25]. In our previous studies, which used the information contained in the spectrograms, similarities were observed in the characteristics of the spectrograms generated by snow avalanches and those associated with other types of GMM (landslides [22], debris flows [23], and lahars [24]). This methodology developed for the study of snow avalanches can be used to detect and characterize GMMs.

In this work, we develop a template with different values of $\beta$ the exponential curves $F\left(t\right)=\mathsf{{\rm K}}{e}^{\beta t}$ with which, once calibrated, the mass movement can be remotely characterized with the information obtained from the spectrogram of a single seismic station. Unlike the amplitude method, which uses more than one station, or the use of DAS, the use of the proposed spectrogram shape method is not affected by the previously mentioned site effects. This is because, being a comparative method in the same location, the site effect that modifies the amplitudes and frequencies remains almost the same.

3. Data

This study uses data obtained at different experimental sites (see the acronym description below) and from three types of GMM: (a) snow avalanches recorded at the VdlS test site (SLF, Switzerland) and Ryggfonn experimental site (NGI, Norway); (b) lahars recorded at the Colima Volcano (CUEIV and UNAM, México); (c) debris flows recorded in the Lattenbach catchment (BOKU, Austria); and (d) the seismic signal corresponding to the Laguna Beach landslide (Caltech/USGS Regional Seismic Network). Different seismic recording stations were used to obtain the seismic data depending on the case.

Table 1. Information on the sites and data acquisition characteristics with references.

Type of GMM	Site	Instruments Sensor	References
Snow avalanches	cavern B VdlS, (SLF, Switzerland)	Mark L-4C-3D (Mark Products) (1 Hz) sensitivity of 277 V/m/s	[26,27,28]
Snow avalanches	Track, Ryggfonn (NGI, Norway)	Mark L-4C-3D (Mark Products) (1 Hz) sensitivity of 277 V/m/s	[29,30]
Lahar	Colima Volcano (CUEIV and UNAM)	Guralp CMG-6TD broadband seismometer (0.03–100 Hz) sensitivity of 2400 V/m/s	[31,32]
Debris flow	Lattenbach catchment (BOKU, Austria)	SM4 (10–180 Hz) sensitivity of 28.8 V/m/s microphone ‘Gefell WME 960H’ (0.5–20 kHz) sensitivity 50 mV/Pa.	[33]
Landslide	Laguna Beach (California, USA)	SCEDC (1–30 Hz) Caltech/USGS Regional Seismic Network	[34,35]

presents the information of the sites and the acquisition instruments and a list of references where detailed information on the GMM used is provided.

For the VdlS test site, we considered the data obtained from cavern B, which was 985 m from the crown (CB1) (Figure 2a). In the case of the Ryggfonn experimental site, we considered the data recorded at the station in the track (TR), which was placed approximately 1450 m from the crown (Figure 3a). In Figure S2, the topographic profiles of the paths of the VdlS and Ryggfonn sites are shown at the same scale for comparison. The data from the Laguna Beach landslide were recorded at station LSS, which is located 23 km from the landslide (Figure 4a). The recording station of the Colima Volcano, México (Figure 5a), was installed ~100 m from the main channel of the Montegrande ravine at ~2100 m.a.s.l. The recording station of the debris flow in the Lattenbach catchment, Austria (Figure 6a), was placed on land beside the channel of the flow at a very short distance (d < 10 m). In addition to the seismic signals generated by the debris flow, infrasound data were also obtained. Data from snow avalanches were obtained at stations consisting of a three-component seismometer, Mark L-4C-3D (Mark Products), with an eigenfrequency of 1 Hz, and a data acquisition system, Reftek–130 (Trimble). The characteristics of the recording stations at the VdlS test site are described in [27,28]. The characteristics of the recording stations at the Ryggfonn experimental site, similar to those of VdlS, are described in [29,30]. The data from the Laguna Beach landslide were recorded at station LSS. This is a 3-component broadband station (1–30 Hz) SCEDC Caltech/USGS Regional Seismic Network [35]. Data from the lahars of the Colima Volcano, México, were recorded via a Guralp CMG-6TD broadband seismometer (0.03–100 Hz frequency range), as described in [31,32]. The complete characteristics of the recording station of the debris flow placed in the Lattenbach catchment (Austria) are described in [33]. The geophone was an SM4 with a frequency range of 10–180 Hz and a sensitivity of 28.8 V/m/s. In addition to the seismic signals generated by the debris flow, infrasound data were also obtained with a nearby microphone, ‘Gefell WME 960H’, with a 0.5–20 kHz frequency range and 50 mV Pa⁻¹. The two types of measurement were obtained with time synchronized with the aim of comparing the suitability and sensitivity of the two types of instrument for monitoring purposes.

The seismic and infrasound data from the experimental sites were collected at 100/200 s.p.s. (50/100 Hz Nyquist frequency) and transformed into ground velocity (m s⁻¹) or pressure (Pa), respectively, via the conversion factors indicated by the manufacturers. In all the cases, the GMM (snow, lahars, or debris) passed over or laterally close to the seismic sensors that were oriented orthogonal to the vertical direction and in the N–S and E–W directions. The seismic signals were bandpass filtered 1–40 Hz with a 4th-order Butterworth filter. Avalanche seismic data were analyzed considering only one seismogram, referred to as the final seismogram (FS). To construct this single-vector seismogram from the three components (Z, N–S, and E–W) that maintain all the seismogram characteristics, the vector properties in 3D space were considered. The construction process is explained in [24]. Only the vector amplitude was considered in the spectrogram. Seismograms of the vertical (Z) component were used for lahars and debris flows. Data from the landslide were recorded at 20 s.p.s. (10 Hz Nyquist frequency) and the E-W component (channel BHE) was used. In the calculation, a short-time fast Fourier transform (FFT) with a Hanning window (length of 0.64 s/1.28 s) and 50% overlap (0.32 s/0.64 s) was used according to the sampling rate. The length and overlap of the time windows for the spectrograms were fixed according to the time series characteristics. Depending on the sampling rate, the resolution in terms of frequency is 0.0156 Hz/0.0078 Hz, and the resolution in time is 0.32 s/0.64 s.

3.1. Snow Avalanches (VdlS and Ryggfonn)

Two sites are considered for this type of mass movement, the VdlS experimental site (Figure 2a and Figure S2) and the Ryggfonn experimental site (Figure 3a and Figure S2). Previously, in [24], the β and K parameters of the SON_spectrogram of different types and sizes of avalanches recorded at site B at VdlS were analyzed (Figure 2a). In conclusion, avalanches can be grouped by similar β and K values according to their type: POW (powder snow), TRANS (transitional), and/or WET (wet snow). A detailed description of each type of avalanche is provided, e.g., [31] and references therein.

Table 2. Information on the signals of snow avalanches recorded at B at the VdlS experimental site. TRANS: transitional snow avalanche. POW: Powder snow avalanche. WET: wet snow avalanche [17]. Sig. Duration (s): Duration of the total signal in seconds. SON Duration (s): Duration of the SON section in seconds. SON/SIG: ratio between the values of the previous columns. Κ (Hz) and β (s⁻¹) with standard error computed values. Last row: averaged values of the different rows.

Event	Sig. Duration (s)	SON Duration (s)	SON/SIG	K (Hz)	β ± σβ (s−1)
TRANS	180	25	0.14	3.10	0.10 ± 0.07
POW	140	30	0.21	1.09	0.14 ± 0.12
WET	60	30	0.50	3.32	0.05 ± 0.02
Average	127 ± 86	28 ± 4	0.28 ± 0.27	2.5 ± 1.7	0.10 ± 0.06

shows the durations of the signals for each avalanche type. The durations of the signals do not exceed 3 min, and the SON_spectrograms range between 14 and 50% of the total duration of the signals. Figure 2 shows, as an example, the time series and the SON_spectrogram of the POW avalanche with the corresponding calculated curves with the values of Table 2 applied to Equation (7). In addition to the durations of the avalanche signals, the β and Κ values are shown in Table 2.

Four avalanches recorded at the Ryggfonn experimental site (Norway, NGI) (Figure 3a and Figure S2) are also included in this study. Two avalanches of each type (dry-dense (d/d) and dry-mixed (d/m)) (NGI classification) are considered (

Table 3. Information on the signals of the snow avalanches at the Ryggfonn experimental site. Event: name of the event. d/m: dry mixed avalanche. d/d: dry dense avalanche. Sig. Duration (s): duration of the total signal in seconds. SON Duration (s): duration of the SON section in seconds. SON/SIG: ratio between the values of the previous columns. Κ (Hz) and β (s⁻¹) with standard error: computed values. Last row: average values of the different rows.

Event	Sig. Duration (s)	SON Duration (s)	SON/SIG	K (Hz)	β ± σβ (s−1)
2004059 d/m	80	40	0.5	0.37	0.91 ± 0.05
2005106 d/d	120	35	0.3	1	0.23 ± 0.11
2007081 d/d	150	30	0.2	0.78	0.14 ± 0.20
2008113 d/m	95	45	0.5	0.04	0.08 ± 0.04
Average	111 ± 53	38 ± 11	0.37 ± 0.25	0.55 ± 0.74	0.13 ± 0.12

). The characteristics of the site and avalanches can be found in [29,30,36] and references therein. The duration of the signals is in the order of 2 min, with d/d being shorter, and the SON section corresponds to 37% of the total duration. In Table 3, the obtained β and Κ values are shown. Considering the duration of the SON_spectrogram sections and the speed of the avalanche front (10 m s⁻¹ for d/d avalanches and 25 m s⁻¹ for d/m avalanches [30]), the first signals of the avalanche come from a distance of approx. 1000 m distance for d/m avalanches and from a distance of 300 m for d/d avalanches (Figure 3a). Thus, the signals of m/d avalanches begin to be recorded at distances from the avalanche farther from the sensor than those of d/d avalanches. These results agree with the length of the avalanche path (1608 m) and with the notion that higher frequencies attenuate with distance more rapidly than lower frequencies do. In the case of the 2,004,059 d/m artificially triggered avalanche, seismic signals are observed in the spectrogram at t = 10 s (Figure 3c), which corresponds to the dashed line in Figure 3a. Since the SBO section is defined as the moment when the avalanche passes over the sensor, we can assume that the avalanche body is over the sensor at 40 s, taking 30 s to reach the sensor from where the avalanche generates enough energy to be recorded.

3.2. Landslide (Laguna Beach)

A significant landslide occurred at the Bluebird Canyon near the town of Laguna Beach in southern California on 1 June 2005, just before 7 a.m. [37,38] (Figure 4a). This landslide occurred in the coverage area of the Caltech/USGS Regional Seismic Network [35]. In [34], the seismic signals of this landslide, recorded at the closest stations, SDD and LLS, at 9 km and 23 km, respectively, from the landslide, were analyzed. The recorded seismic signals with a duration of approx. 11 h show 4 blocks of energetic arrivals of different durations, which are better observed at station LSS-BHE, and interpreted as being generated by the landslide (Figure 5 in [34]). This indicates that the landslide occurred in at least in 4 stages. Figure 4b,c show the seismic signal of the first 11,000 s (3 h) and its spectrogram. We assume that the first section of 2000 s, indicated between the dashed lines in Figure 4b,c, with a smaller amplitude than the later sections, corresponds to the first stage, the beginning of the landslide. Figure 4d,e show the seismic signal and spectrogram of this part, with the calculated exponential curve corresponding to the Κ and β values (Equation (7)). The SON_spectrogram section has a duration of 800 s and corresponds to 32% of the total duration of the signal. The obtained Κ and β values combined with the durations of the signal, SON section, and its ratio are shown in

Table 4. Information on the seismic signals of the Laguna Beach landslide. Signals: type of measurement. Sig. Duration (s): Duration of the total signal in seconds. SON Duration (s): Duration of the SON section in seconds. SON/SIG: ratio between the values of the previous columns. Κ (Hz) and β (s⁻¹) with standard error: computed values.

Signals	Sig. Duration (s)	SON Duration (s)	SON/SIG	K (Hz)	β ± σβ (s−1)
First part	2500	800	0.032	0.49	0.00391 ± 0.00001

.

3.3. Lahars

The lahars of the Colima Volcano were monitored by different instruments installed ~100 m away from the Montegrande ravine (Figure 5). The site and the data acquisition system are described in [28], and we considered the same data used in this study. In that work, the evolution of the frequency content of the seismic signals was used to discern the different seismic sections connected with the different phases of the lahars (front, body, and tail). Spectrograms and different frequency bands were analyzed for this purpose. Lahars that occurred on 25 June and 15 September 2012 and 11 June 2013 and 24 July 2013 are analyzed in this work. In [39], the activity of the volcano during that period was examined. The characteristics of the 15 September 2012 lahar are described in [40].

Table 5. Information on the signals of the Colima Volcano studied lahars. Event: name of the lahar. Sig. Duration (s): duration of the total signal in seconds. SON Duration (s): duration of the SON section in seconds. SON/SIG: ratio between the values of the previous columns. Κ (Hz) and β (s⁻¹) with standard error computed values. Last row: averaged values of the different rows.

Event	Sig. Duration (s)	SON Duration (s)	SON/SIG	K (Hz)	β ± σβ (s−1)
20120625	3000	700	0.23	1.000	0.0040 ± 0.0020
20120915	4500	1000	0.22	0.188	0.0035 ± 0.0011
20130611	5000	800	0.16	1.749	0.0024 ± 0.0006
20130724	5000	1500	0.30	0.920	0.0030 ± 0.0003
Average	4375 ± 1639	1000 ± 616	0.23 ± 0.01	0.96 ± 1.11	0.003 ± 0.001

shows the seismic signal durations of the lahars, which were approximately half an hour. The larger lahar had a duration of approximately 5000 s (24 July 2013), with an SON section of approximately 1500 s. The SON_spectrograms in the analyzed lahars correspond to 16 to 30% of the total duration of the signal (Table 5). Table 5 also shows the β and K values obtained via the proposed method. Figure 5 shows the seismogram and SON_spectrogram of the first 2000 s of the 2,013,0724 lahars with the F(t) curves (Equation (7)) superimposed and calculated with the β and K values indicated in Table 5.

3.4. Debris Flows

The data were obtained in the Lattenbach Torrent in west Tyrol (Figure 6a), which was monitored by different instruments. In [33], a description of the Lattenbach site and analyzed debris flows was presented. In Figure 6, the time series and SON_spectrograms of the infrasound (b and c) and geophone (d and e) signals of the debris flow with the curves according to

Table 6. Information on the seismic and infrasound signals of the Lattenbach debris flows. Signals: type of measurement. Sig. Duration (s): duration of the total signal in seconds. SON Duration (s): duration of the SON section in seconds. SON/SIG: ratio between the values of the previous columns. Κ (Hz) and β (s − 1) with standard error: computed values.

Signals	Sig. Duration (s)	SON Duration (s)	SON/SIG	K (Hz)	β ± σβ (s−1)
seismic	2000	78	0.04	1.01	0.01 ± 0.01
Infrasound	3000	200	0.07	0.04	0.017 ± 0.005

and Equation (7) are shown. Because the seismic data were recorded in a 10 Hz vertical geophone, unfortunately, no frequency below this value can be obtained. However, the infrasound data have a lower frequency content, which allows us to better detect the increasing shape of the SON sections. In Table 6, the durations of the signals recorded by the geophone and infrasound sensor are shown. The duration of the infrasound signals was in the order of one hour, and the SON_spectrogram corresponded to 7% of the total duration. The corresponding values for the seismic data are not considered because of their inadequate quality. Table 6 also indicates the values of β and K obtained from both types of data. The infrasound values were used to calculate the curves plotted on the SON_spectrograms (Figure 6). Note that the curves are also valid for the seismic data. This result confirms the existence of coupling between infrasound and seismic data [41,42]. The comparison between the time series and the spectrograms of the signals from both sensors (speed of the ground and variations in pressure) were carried out in [27], indicating an early detection of the infrasound with respect to the seismic data. However, considering our previous comparison of the SON_spectrograms, this result could be a consequence of the cutoff of the low frequencies of the geophone due to the high natural frequency of the sensor. The duration of the SON_spectrogram indicates that the first signals of the debris flows were generated approximately 1400 m from the sensor (Figure 6a). To obtain this value, we considered that the speed of infrasound transmission in the air (approx. 340 m s⁻¹), although not instantaneous, is two orders of magnitude greater than the speed of the event (7 m s⁻¹ [27]). This value agrees with the topographic characteristics of the area.

4. Methodology: Theoretical Assumptions

This section is devoted to the theoretical justification of why we can identify the exponential shape in the SON_spectrogram section and link it to the seismic waves propagation properties in the ground. In [24], a method to estimate the parameters that define the exponential shape of the SON_spectrogram was presented (Figure 1b). GMMs are moving sources of seismic energy that generate a seismic wave field. The propagation of seismic energy in the earth is disturbed by different factors, such as seismic wave attenuation [43], which affects the amplitude and frequency content of the generated wave field and must be considered when its records (e.g., seismograms) are analyzed.

As stated above, one distinctive feature depicted in the SON_spectrograms is the increase in the amplitude of the high frequencies in an exponential shape over time. We assume that this increase is due to wave propagation in a medium, which is due mainly to the attenuation effect (anelastic/intrinsic and geometrical spreading, e.g., [44,45]) of the seismic waves, and possibly increases in mass by entrainment as the mass approaches the sensor. This increasing pattern could be attributed to the Doppler effect.

However, the quantification undertaken in [23], in the case of snow avalanches, rejected this because of the difference of two orders of magnitude between the speed of the avalanche (∼10 m s⁻¹) and that of the seismic waves traveling across the ground (∼10³ m s⁻¹). This can also be valid in the case of debris and lahars given the similar order of their velocity values (∼10 m s⁻¹) (e.g., [28,46]).

Synthetic Spectrogram

The flowchart in Figure 7 shows the procedure of interpretation of the shape of the SON_spectrograms following the theoretical justification. To reproduce this shape, a theoretical spectrogram was created via Equation (1) of the theoretical decay of the seismic amplitude and frequency content with distance. Equation (1) has already been used by other authors for localization via the amplitudes of seismic signals (e.g., [4,47,48]). This equation considers both attenuation by geometrical spreading and intrinsic attenuation as they move away from the energy source. The source of energy is assumed to be the descending mass. The equation estimates the wave amplitude, $\mathbf{A}\left(\mathit{t},\mathit{f}\right)$ given at a point, at r′ distance from a source at r′ = 0, with an amplitude, $\mathit{A}\mathit{S}(f_s,t_s)$, given as a function of $t_s$ (time) and $f_s$ (frequency), as follows:

$$\mathbf{A}\left(\mathit{t},\mathit{f}\right)=\mathit{A}\mathit{S}(f_s,t_s){\displaystyle \frac{1}{{\left(2\mathsf{\pi }\mathit{h}{\mathit{r}}^{\prime }\left(t\right)\right)}^{1/2}}}{\mathit{e}}^{-\mathit{\alpha }\left(\mathit{f}\right){\mathit{r}}^{\prime }\left(\mathit{t}\right)},$$

(1)

where $\mathit{A}\mathit{S}(f_s,t_s)$ is the amplitude that a hypothetical sensor would measure at r’ = 0; $\mathit{A}\mathit{S}(f_s,t_s)$ can be considered to last a time interval. The quotient term corresponds to the geometrical spreading. This is a decreasing function as distance r’ increases from the source. This equation considers the energy flux of the seismic wave through a cylindrical surface of height h and radius r’.

The exponential term corresponds to the intrinsic attenuation, where:

$$\mathit{\alpha }\left(\mathit{f}\right)={\displaystyle \frac{\mathsf{\pi }\mathit{f}}{\mathit{Q}\left(\mathit{f}\right)\mathit{c}\left(\mathit{f}\right)}}$$

(2)

and $h=\mathrm{k}\mathit{c}\left(\mathit{f}\right)/\mathit{f}$. The factor k is associated with the skin depth [48,49]. $\mathit{Q}\left(\mathit{f}\right)$ ground quality factor (dimensionless) and $\mathit{c}\left(\mathit{f}\right)$ is the phase velocity of the surface waves. The intrinsic attenuation affects the frequency content of the wave.

Since our attention is focused on the signals recorded at the sensor, at r′ distance from a source, Equation (1) can also be interpreted as the variation in the frequency content at the sensor as the source approaches it. Then, we can consider the source approaching the sensor according to r(t) (with r = 0 at the sensor). This implies a change in the coordinate system, with R being the distance between the source and the sensor at a fixed time, $t$₀. A scheme of the dual-distance coordinate system is shown in Figure S4. Since $\mathit{R}=\left|\mathit{r}\left(t\right)\right|+\left|{\mathit{r}}^{\prime }\left(t\right)\right|$, ${\mathit{r}}^{\prime }\left(t\right)$ can be expressed as $\left|{\mathit{r}}^{\prime }\left(t\right)\right|=\mathit{R}-\left|\mathit{r}\left(t\right)\right|$, where ${\mathit{r}}^{\prime }\left(t\right)$ and $\mathit{r}\left(t\right)$ are the complementary distances along the path. Note that, as the point source approaches the sensor, r becomes 0 and r’ = R. Moreover, r increases in time because −R < r < 0, and Equation (1) becomes the following increasing function:

$$\mathbf{A}\left(\mathit{t},\mathit{f}\right)=\mathit{A}\mathit{S}(f_s,t_s){\displaystyle \frac{1}{{\left(2\mathsf{\pi }\mathit{h}\left(\mathit{R}-\left|\mathit{r}\left(t\right)\right|\right)\right)}^{1/2}}}{\mathit{e}}^{-\mathit{\alpha }\left(\mathit{f}\right)\mathit{R}}{\mathit{e}}^{\mathit{\alpha }\left(\mathit{f}\right)\mathit{r}\left(\mathit{t}\right)}$$

(3)

which has its maximum at r = 0.

Assuming that ${\mathit{v}}_{\mathit{m}}$ (t) is the mass speed, r(t) can be expressed as r(t) = ${\mathit{v}}_{\mathit{m}}$ (t) t, and Equation (3) can be expressed as

$$\mathbf{A}\left(\mathit{t},\mathit{f}\right)=\mathrm{K}\left(\mathit{t},\mathit{f}\right){\mathit{e}}^{\mathit{\beta }\mathit{t}},$$

(4)

where:

$$\mathbf{A}\left(\mathit{t},\mathit{f}\right)=\mathit{A}\mathit{S}(f_s,t_s){\displaystyle \frac{1}{{\left(2\mathsf{\pi }\mathit{h}\left(\mathit{R}-\left|\mathit{r}\left(t\right)\right|\right)\right)}^{1/2}}}{\mathit{e}}^{-\mathit{\alpha }\left(\mathit{f}\right)\mathit{R}}{\mathit{e}}^{\mathit{\alpha }\left(\mathit{f}\right)\mathit{r}\left(\mathit{t}\right)}$$

(5)

and

$$\beta =\mathit{\alpha }\left(\mathit{f}\right){\mathit{v}}_{\mathit{m}}\left(t\right)$$

(6)

In these expressions, we use β instead of β(t, f) for simplicity. Note that the equation indicates the dependence of β on the intrinsic attenuation, $\mathit{\alpha }$, and on the GM speed, ${\mathit{v}}_{\mathit{m}}\left(t\right)$. If β and $\mathit{\alpha }$ are known, an estimate of the mass velocity can be obtained.

Equation (4) represents the variation of the amplitude in time and frequency (or the variation in time of the frequency content) of the signal and is an increasing function in time, t. This increase depends on $\mathrm{K}\left(t,f\right)$, which depends on the decreasing function on time, |r(t)|^1/2, and on the exponential function of t modulated by β, which depends on the ground characteristics, which are frequency dependent. The factor $e^{-\mathit{\alpha }\left(\mathit{f}\right)\mathit{R}}$ is a constant value of r and t.

Figure 7a shows the synthetic spectrogram of a theoretical SON_spectrogram section via Equation (5). In this case, a unitary impulse of amplitude $(\mathit{A}\mathit{S}\left(f_s,t_s\right)$ = 1 u) is considered. The terrain values used, $\mathit{Q}\left(\mathit{f}\right)$ and $\mathit{c}\left(\mathit{f}\right)$, are those obtained in [32]. The velocity of the approaching mass to the sensor is assumed to be ${\mathit{v}}_{\mathit{m}}$ = 10 m s⁻¹ and the factor k = 0.25. Notice the similarity in the shape of this theoretical spectrogram with those of the SON_spectrograms obtained from the seismic signals (e.g., Figure 1b). The time-increasing shape of the amplitudes in the SON_spectrogram is determined by the dependence on the frequency of the β exponent coefficient (Equation (6)). This would have the same effect as an increase in the GM. Since $\mathit{\alpha }\left(\mathit{f}\right)$ is independent of time and distance, its role in Equation (6) is a constant of multiplication over time. Having established the relationship between a synthetic spectrogram and one generated from a seismogram, we are now interested in estimating the K and β values from a spectrogram of measured data. These values of the SON_spectrogram allow us to obtain information on the physical process of mass movement since they are related mathematically to it (Equations (4)–(6)). In [24], a sequential procedure was used to obtain these parameters from the spectrograms, and the function

$$\mathrm{F}\left(t,f\right)=\mathsf{{\rm K}}{e}^{\mathsf{\beta }t}$$

(7)

is presented.

5. Results

In

Table 7. Summary of the averaged values of the information of the studied signals. EVENTS: event type. SD ± σSD (s): average duration of the total signal in seconds. SOND ± σSOND (s): averaged of SON_spectrogram section duration in seconds: SON/SIG ± σ: averaged ratio between the values of the previous columns and the average values of Κ and β with standard error.

EVENTS	SD ± σSD (s)	SOND ± σSOND (s)	SON/SIG ± σ	K ± σK (Hz)	β ± σβ (s−1)
Snow av. (VdlS)	127 ± 86	28 ± 4	0.28 ± 0.27	2.50 ± 1.74	0.10 ± 0.07
Snow av. (RYGGFONN)	111 ± 53	38 ± 11	0.37 ± 0.25	0.55 ± 0.74	0.13 ± 0.12
Landslide (Laguna Beach)	2000	800	0.4	0.49	(3.91 ± 0.01) 10−3
Lahars	4375 ± 1639	1000 ± 616	0.23 ± 0.01	0.96 ± 1.1	0.003 ± 0.001
Debris (I)	2000	200	0.10	0.039	0.017

, the average values of the parameters obtained for each type of mass movement are presented. Differences among them are evident. The duration of the avalanche signals (no more than 3 min) is much shorter than that of the lahar, landslide, and debris flow signals, with durations in the order of hours and half an hour, respectively. The duration of the SON_ spectrogram section (SOND) of the avalanches is smaller than those of the debris, lahars, and landslides, by one and two orders of magnitude, respectively. The SON/SIG ratio of the debris flow is smaller than those of the other GMMs, although we are aware that we have only data for one event. With respect to K and β, independent of the type of mass movement, these two parameters can be used to characterize the increasing shape of the SON_spectrogram section. The increasing shape connected to a mass movement approaching a sensor is described by these parameters according to the type of event (Table 7). The average β value for lahars is approximately 0.003 s⁻¹, whereas it is approximately 0.0004 s ⁻¹ for the landslide, which is two orders of magnitude lower than that for avalanches (approximately 0.1 s⁻¹). Those of the debris flows are approximately 0.01 s⁻¹. The β value, the coefficient of the exponential (Equation (7)), is responsible for the increasing shape. The effect of K is a shift of the curve along the axes, as can be seen in Figure S5. With respect to the average values of β for avalanches, regardless of the site (VdlS or Ryggfonn, Figure S2), these values are similar.

For the Laguna Beach landslide, the spindle shape observed in the time series and the shape of the spectrogram are similar to those of the other GMMs previously analyzed (e.g., [24]). Although the situation is different from those proposed in the theoretical assumption presented here, where the increasing shape is essentially justified by the effect of the intrinsic attenuation and the geometric spreading without considering mass increase, a possible explanation for this similar behavior would be that, in the case of a distant GMM, the increase in mass (amplitude in source $(\mathit{A}\mathit{S}\left(f_s,t_s\right),\mathrm{E}\mathrm{q}\mathrm{u}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}(5$)) would be responsible for the increasing shape, which is dominant over the attenuation effect. The geometrical spreading over large distances cannot affect the amplitudes since the size of the landslide is negligible in relation to the distance to the sensor (Figure 4a), and at the distances considered, the high frequencies have already been attenuated, with only the low frequencies remaining. In Figure 8a, the curves obtained with the average parameter values of the events investigated are presented simultaneously. To better interpret these results, Figure 8b shows the same curves but with time zero corresponding to when the mass movement is over the sensor. The curves indicate that the appearance of high-frequency seismic energy with time is faster in avalanches than in debris flows and slowest in lahars and landslide. In the case of mass movement, the interval for its detection in a nearby sensor will be longer for lahars and landslides than for the other GMMs, of the order of 1200 s in the former and 200 s or 20 s in the latter (Figure 8 and Table 7). In other words, there will be more time for the alert if a lahar or landslide occurs than in the case of a snow avalanche, which occurs at a faster rate. Note that, in the case of the landslide, the curve does not reach the sensor and stops at −500 s, which is a consequence of the cutoff of the high frequencies due to sampling (20 s.p.s.). We have shown that, in the case of the GMM analyzed, the different characteristics of each of them are reflected in the shape of the spectrogram (K and β parameters). Since the exponential form is present, not only in the GMMs analyzed but also in general, under the established conditions, we believe that this behavior and methodology may be useful for seismic characterization in other cases.

We have therefore created a template of curves with different K and β parameters that can detect and characterize the type of events at a distance to aid professionals involved in risk management. However, to simplify the process, the role of these parameters in this template is explained. Previously, we mentioned that the effect of K is a shift along the axes. Figure S4 shows the curves created using the values of K and β for snow avalanches at Ryggfonn (Table 3) and those created using the same β value and setting K = 1 Hz, where it is shown that the effect of K is a shift only along the time axis (horizontal black arrows). Considering the effect of K, a template can be created for different values of β (K = 1 Hz). Figure 9 shows a template of curves created with different values of β (K = 1 Hz) on the SON_spectrogram of avalanche POW (Figure 2b) recorded by sensor B at VdlS (Figure 2a). Note that the greater β is, the greater the increase in slope, or in other words, the faster the appearance of energy at high frequencies in the SON_spectrogram. To obtain the β value of the observed spectrogram via the template, it is not necessary to develop the algorithm or the procedure presented previously [24]. Once the template is calibrated with the range of β values according to the type of mass movement and the specific site, a comparison of the spectrogram with the template at the same scale can help in the classification within each type of mass movement. Figure 9 shows the comparison of the template with the SON_spectrogram section. The curve with a value of β = 0.14 s⁻¹ (the dotted line over the white line) fits the shape of the SON_spectrogram. This is just the value we obtained for the POW avalanche (Table 2). Figure 10 shows the use of the template for the SON_spectrogram of the 2,012,0915 lahar at the Colima Volcano as an example to obtain the β coefficient. Once the template was constructed with the values corresponding to the GMM, it slides over the spectrogram until the curve that fits the curvature of the SON_spectrogram is obtained. The value obtained for β is in the range of 0.0035–0.004 s⁻¹, similar to that obtained from the procedure (Table 5).

6. Conclusions

The different GMMs (snow avalanches, landslides, lahars, and debris flows) display similar behavior in the spectrograms in terms of the shape of their initial section (SON_spectrogram) which corresponds to the beginning of the gravitational mass movement when it approaches a seismic sensor. The proposed analysis methodology allows us to obtain the propagation characteristics of different types of GMM by analyzing the spectrograms of the generated seismic signals. We present a template that helps professionals characterize the exponential shape of the SON_spectrogram of GMMs with a single seismic record (professional seismometer) whenever the moving mass approaches the sensor, is close to it, or passes over it. The template must be calibrated at each site considering the site effects (or site response), terrain characteristics, and the type of mass movement. This template is the result of the link between the propagation properties of seismic waves and the fitting of an exponential function to the shape of the SON_spectrogram. The exponent parameter is responsible for the increasing shape.

In addition, we obtained characteristic values of the exponential parameter for each type of GMM. The average value is approximately 0.003 s⁻¹ for lahars, 0.004 s⁻¹ for landslides, and one order of magnitude greater (0.017 s⁻¹) and one order of magnitude lower for debris flows than that of the avalanches (0.12 s⁻¹). Comparing the SON_spectrogram of an event with the appropriate template (parameter value) allows us (just by superimposing it and sliding it graphically) to classify the event. In addition, once the parameter value is determined, it is possible to establish the characteristics of the GMM as its speed increases [24]. However, to determine this approach speed, it is necessary to know the seismic properties of the terrain. For a successful application, the frequency band of the recorded signal must be wide enough (e.g., 2–40 Hz) to be able to construct the exponential curve. In the case of large distances, although the signal can be detected, the application of the template is more problematic.