# On the use of earthquake multiplets to study fractures and the temporal evolution of an active volcano

G. Poupinet, A. Ratdomopurbo, O. Coutant

► To cite this version:

G. Poupinet, A. Ratdomopurbo, O. Coutant. On the use of earthquake multiplets to study fractures and the temporal evolution of an active volcano. Annals of Geophysics, 1996, 39 (2), 10.4401/ag-3968. insu-03448565

HAL Id: insu-03448565

https://insu.hal.science/insu-03448565

Submitted on 25 Nov 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# On the use of earthquake multiplets to study fractures and the temporal evolution of an active volcano

Georges Poupinet(¹), Antonius Ratdomopurbo(¹)(²) and Olivier Coutant(¹)

(¹) LGIT-IRIGM, Université Joseph Fourier and CNRS, Grenoble, France

(²) Volcanological Survey of Indonesia, MVO, Jl Cendana, Jogyakarta, Indonesia

## Abstract

Multiplets, i.e. events with similar waveforms, are common features on active volcanoes. The seismograms of multiplets are analyzed by cross-spectrum techniques: this procedure improves by a factor of about 10 the precision of differential P-arrival times and therefore the accuracy of the relative location of earthquakes. Long period events which cannot be located because of the impossibility to pick up P-waves on individual seismograms can be located with a precision of about 10 m. Such a precision permits fault planes to be mapped inside a volcanic edifice and the azimuth and strike of fractures to be defined. Seismograms of the two events (of a doublet) that occur on different dates are analyzed by the Cross Spectrum Moving Window technique (CSMW) for measuring the time delay between waves in the coda. The pattern of the delays in the coda is a function of the temporal changes of seismic velocity that occurred inside the volcano during the time interval that separates the two events of a doublet. We illustrate the potential of the doublet technique for detecting temporal changes inside a volcano by performing computations of synthetic seismograms. The case of a dyke injected inside the volcano is considered as well as that of the replenishment of a superficial magma chamber and of a general increase in velocity in the summit of the volcano. Data from Merapi volcano (Indonesia) illustrate a possible temporal velocity change inside the volcano several months before the 1992 eruption.

**Key words** volcanic earthquakes – relocation of events – temporal monitoring

## 1. Introduction

Among tectonic earthquakes, similar waveform events are sometimes found. A pair of earthquakes with similar waveforms is called a doublet. Multiplets (families of similar earthquakes) are frequently observed on volcanoes (Okada et al., 1981; Frémont and Malone, 1987; Got et al., 1994). The trace similarity of two seismograms of a doublet means that their hypocenters are spatially near, that their sources are equivalent and that the ray trajectories from sources to receivers are about the same (within much less than one wavelength).

The spatial gap between sources is smaller than the source-receiver distance: a few meters compared to kilometers. Identification of seismic multiplets can be performed by visual observation of the seismograms, finding trace similarity. Another way to find multiplets is to build a database with all the events recorded (for instance a short duration record for each source-station pair) and to compute all the cross-correlations between these events. From a cross-correlation matrix the events that have the same arrival time differences are extracted (Got et al., 1994). Figure 1 is an example of the events of a multiplet recorded in station PUSV on Merapi volcano (Indonesia).

Several studies on volcanic doublets have improved the precision in the relative location of events. Doublet analysis is also useful for
<figure><img src="image_2.png" /><figcaption>Fig. 1. Example of several shallow volcano-tectonic earthquakes (VTB) recorded in station PUSV on the Merapi volcano. These events are called a multiplet.</figcaption></figure>

measuring temporal changes in seismic velocity, for instance to detect a velocity change in the vicinity of a fault that may precede or follow a large magnitude earthquake due to stress redistribution: velocity changes of a few tenths of a percent can be detected (Poupinet et al., 1984).

## 2. Cross-spectrum analysis of doublets

### 2.1. Cross-spectrum analysis of P-waves and S-waves

Any location of earthquakes requires the precise timing of at least P-arrival times. Using doublets, we compare the waveforms of the P-waves of different earthquakes in each station and we measure the arrival time differences

between the two events by a standard cross-correlation technique. First, the seismograms are aligned with a precision of one sample, usually 10 or 5 ms. A window is applied on the P-wave train: usually 128 points (1.28 s) tapered by a Hanning function are adequate. In a given window, the phase spectrum of the cross-correlation between two events, $\phi$, is directly proportional to the frequency $f$ and to the time delay $\Delta\tau$ between the two signals

$$ \phi(f) = 2\pi\Delta\tau f. $$

Thus the delay $\Delta\tau$ between the two traces is computed from the slope $b$ of the phase of the cross-spectrum

$$ \Delta\tau = b/(2\pi). $$
The phase error $\sigma$ is estimated

$$ \sigma^2 = \frac{1}{2B_w T} \left( \frac{1}{|C|^2} - 1 \right) $$

where C is the coherency level, $B_w$ and T are the bandpass and the duration of the signal respectively. The slope b is estimated by a weighted linear regression that takes into account the phase error $\sigma$ as a weighting factor. Figure 2 shows the basic processing of a dou-

blet on data recorded on Merapi volcano by the Merapi Volcano Observatory (MVO) in Jogyakarta: the seismograms of the two events constituting a doublet are plotted on top of fig. 2. The coherency and the phase spectrum of the cross-correlation of the P-wave train (box) are plotted as a function of frequency in the bottom of fig. 2. These events were recorded with a Mark Product L4C-1 Hz seismometer; the seismic signal is telemetered to MVO where it is digitized and recorded on a PC microcomputer at a rate of 100 samples per second with IASPEI software (Lee, 1989). It

<figure><img src="image_3.png" /><figcaption>Fig. 2. Cross-spectral analysis of the P-wave trains (box) of two VTB earthquakes recorded in station PUSV. The coherency of the boxed signal is good until 16 Hz. The phase spectrum of the cross-correlation is linear: the slope of the phase spectrum gives the delay between the two P-wave trains. The digitization rate is 100 points per second but a precision of the order of 0.002 s is achieved in the relative timing of these two P-waves. A precision of up to 10 times smaller than the digitization rate can be achieved.</figcaption></figure>
should be noted that this processing can give a relative timing precision of about one tenth the sampling rate because it corresponds to an adjustment of the phases of the signals.

## 2.2. Relocation of doublets

Seismic events are located in routine analysis because it is possible to pick up their P-waves and sometimes their S-waves. Several authors have relocated doublets using either a

cross-spectral technique in the frequency domain as described earlier or a cross correlation technique in the time domain (Poupinet et al., 1982, 1984; Scherbaum and Wendler, 1986; Frémont and Malone, 1987; Deichman and Garcia-Fernandez, 1992; Got et al., 1994). Both techniques are equivalent. The relocation is performed by minimizing

$$ \Sigma(\Delta\tau^i_{\text{observed}} - \Delta\tau^i_{\text{computed}})(\Delta X, \Delta Y, \Delta Z, \Delta h)^2 \\ i = 1, n $$

<figure><img src="image_4.png" /><figcaption>**Fig. 3a.** Hypocenters of the VTB multiplet presented in fig. 1. The top plot is an epicenter map and the bottom one is a vertical cross-section. The origin of the coordinates is in the crater of Merapi at an elevation of 2925 m. These events are located beneath the south flank of the volcano at a depth of less than 1 km from the summit.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 3b.** Relocation of the VTB multiplet presented in fig. 1 using the doublet technique. The top plot is an epicenter map and the bottom one is a vertical cross-section. The precise relocations define a fracture near the summit of the volcano striking E-W and nearly vertical. The error on location is of the order of 10 m. Note the difference in scale between fig. 3a and 3b.</figcaption></figure>
where *i* is the index of stations. $\Delta\tau^i$ is the observed difference in arrival times between the two *P*-waves for station *i*. $\Delta T^i$ is the difference between the travel times of *P*-waves computed for two events whose hypocenters differ in coordinates by $\Delta X$, $\Delta Y$, $\Delta Z$ and in origin time by $\Delta h$. This function is dependent on the velocity in the vicinity of the hypocenters and on the azimuth and incidence angles of the stations. The minimum of this function can be found either by exploring the 4D space around an initial hypocenter or by a least-square inversion (see Got *et al.*, 1994 for example).

An illustration of the relocation of VTB type earthquakes from Merapi is shown in fig. 3a,b: 8 events were originally located within 1 km (fig. 3a) by the MVO seismic array which is composed of 6 stations. The same

events are concentrated within a few dozen meters after relocation by the cross-spectral technique (fig. 3b). The precision of the relative location is of the order of 10 m. After relocation these events are aligned and define a fault segment or a crack oriented in the east-west direction. The direction of this crack can then be compared to deformation measurements.

## 2.3. *Cross-Spectrum Moving-Window technique and the analysis of the entire seismogram*

The cross-spectral determination of the time delay between the two seismograms can be performed for any window along the seismo-

<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Cross-Spectrum Moving-Window analysis of an LF doublet (low frequency). The first event occurred on August 30, 1991, and the second on the day after. The coda of the LF event is dominated by several separate phases with different arrival delay times.</figcaption></figure>
gram. Therefore a narrow window is shifted along the entire seismograms. The time delay is plotted as a function of the running time along the seismogram: we call this technique the Cross-Spectrum Moving-Window (CSMW) technique. The CSMW technique gives an accurate measurement of the time delay between all waves along the seismogram, without knowing their exact nature and ray path.

Figure 4 shows the application of the CSMW technique to a low frequency doublet recorded on Merapi. A window of 2.56 s was taken in the analysis. The first arrivals in both events were aligned within one sample. We observe different time delays for individual phases within the coda. Steps in the CSMW diagram help to define the occurrence of waves of different type and ray path.

### 3. Temporal variations of physical properties inside a volcano

A temporal variation of seismic velocity induces travel-time differences for waves traveling through the perturbed medium. If it occurs in the time period between the two events of a doublet, these differences can be measured with the CSMW technique. Any change in arrival time of a wave in the entire seismogram, including a back-scattered wave, gives a travel-time delay $\Delta\tau(t)$. Furthermore, different size and shape of the perturbed volume produce different patterns of $\Delta\tau(t)$ curves. Thus, when CSMW accuracy allows it, we use $\Delta\tau(t)$ measurements to map velocity variations inside the medium. We illustrate these possibilities by showing synthetic data corresponding to four simple but characteristic temporal variations that may occur inside a volcano. We always compute the synthetic seismogram for the unperturbed case then the synthetic seismogram for the perturbed case. We consider these seismograms as the two events of a doublet. We apply the CSMW technique to enhance the differences in waveform between these seismograms and to show what could be measured from surface observations.

#### 3.1. Examples of synthetic CSMW technique plots

##### 3.1.1. A global change in velocity in the entire summit of the volcano

In this case, velocity changes in the entire medium. All waves that propagate within the medium are affected by the velocity change. The monitoring of this velocity variation is based on the intrinsic property of the coda of local seismograms: the coda is essentially the superposition of S back-scattered waves that are radiated by propagating in a heterogeneous medium (Aki and Chouet, 1975). The longer the path traveled by the scattered wave, the later the wave is recorded on the seismogram. Assuming a constant velocity in the medium and a temporal change of velocity from $V_1$ to $V_2$, with $\langle V \rangle = \sqrt{V_1 V_2}$, the travel-time delay $\Delta\tau$ is expressed as:

$$ \Delta\tau = \frac{x}{V_1} - \frac{x}{V_2} = \frac{x\Delta\nu}{V_1V_2} = \frac{x}{\langle V \rangle} \frac{\Delta\nu}{\langle V \rangle} = \langle T \rangle \frac{\Delta\nu}{\langle V \rangle} $$

where the slope of $\Delta\tau$ versus propagation time $\langle T \rangle$ curve is equal to the relative S-velocity variation.

This technique for measuring a temporal change in S-velocity precisely is illustrated in fig. 5. We compute P-SV elastic waves in a heterogeneous random medium by a finite differences method (Virieux, 1986). Propagation through this kind of medium generates scattered waves, the coda waves, as shown on the seismograms. In fig. 5 (bottom), the continuous line represents a vertical seismogram at the surface for the reference (unperturbed) model; the dashed line represents the seismogram for the same source and station when velocity changes by one per cent (perturbed model). We clearly observe a time shift of coda waves. Figure 5 (upper part) is the Cross-Spectrum Moving-Window measurement of the delay between the perturbed seismogram and the reference seismogram. The slope $\Delta\tau/\Delta T$ of the delay as a function of running time is constant along the seismogram. The measure of this slope yields an estimate of the velocity variation of about 1% ± 0.5%.
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Effect of a constant velocity perturbation in the entire medium using a finite difference modeling of the P-SV waves propagation for a random medium with 10% velocity random fluctuations. The vertical seismograms recorded at the surface are plotted before (continuous line) and after the velocity decrease (dashed line) in the bottom of the figure. The size of the window used for the computation of the delay is shown in the lower left corner. We clearly observe a shift of waves in the coda. The upper diagram shows the Cross-Spectrum Moving-Window analysis obtained for the two seismograms below. The CSMW plot is a measurement of the time delay between the waves along the entire seismogram. The velocity perturbation between them is 1 per cent. In this example the time delay is proportional to the velocity change of S-waves because the coda is essentially back-scattered S-energy.</figcaption></figure>

### 3.1.2. A localized velocity variation due to a horizontal dyke intrusion

Let us consider that a horizontal dyke is injected at depth in the time interval between the two events of a doublet. This is represented by a horizontal layer where velocity has changed between doublets. The perturbed medium is not large enough to trap scattered waves so we

cannot estimate the velocity perturbation from coda wave analysis. However, the pattern observed on CSMW plots (fig. 6a,b) allows us to locate this variation. All waves passing through the slow layer (the injected dyke), such as scattered or reflected waves generated below it, are delayed and exhibit a delay with respect to waves propagating in the unperturbed medium. By picking up the arrival time of the $\Delta\tau$ dis-
<figure><img src="image_8.png" /></figure>

**Fig. 6a,b.** Effect of the intrusion of a horizontal dyke. A dyke intrudes a horizontally layered model: it is represented by a decrease in velocity in a layer at a depth of 100 m (a). The corresponding seismograms (before and after the injection of the dyke) and the CSMW diagram which measures the arrival time differences between the various waves of the seismograms are plotted in (b). The two seismograms are so similar that it is extremely difficult to visualize their differences except by using the CSMW technique.

continuity recorded in different stations, we can migrate this $\Delta\tau$ discontinuity to its location assuming a constant and known velocity. We illustrate this case by computing acoustic seismograms with a reflectivity method for horizontally layered models. One of the layers has a 2% velocity variation (from 4.8 to 4.7 km/s) and the source is located at a depth of 75 m. The $\Delta\tau$ discontinuity is picked up at 31 ms (27 ms + window analysis length) on a surface receiver. This travel-time corresponds roughly to a path length of 135 m which in turn, gives us a 100 m depth for the horizontal dyke location.

### 3.1.3. A velocity perturbation inside a superficial magma chamber

We introduce now a velocity change affecting only a limited volume sampled by seismic waves: some waves travel through the perturbed medium and others do not. This can of

course create very complex patterns in a CSMW plot. A simple case is that of direct waves propagating in the unperturbed medium while coda waves with longer paths have been scattered in the perturbed medium. The delay $\Delta\tau$ between the reference and the perturbed seismogram remains zero for direct *P*- and direct *S*-wave trains. The delay increases as a function of time in their respective codas (*P*-coda and *S*-coda). This case is illustrated in fig. 7a,b. The velocity has changed in a limited shallow volume similar to a magmatic chamber denoted by the shaded gray area. Earthquake doublets have been recorded at the surface before and after the velocity changes. The doublets are located above the magma chamber. This simulation has been carried out using a reflectivity code and the coda is generated by multiple reflections in the layers. The CSMW plot shows that the direct *P*- and *S*-wave propagation times are unchanged while velocity variations are recovered from the respective coda delays $\Delta\tau(t)$.
<figure><img src="image_9.png" /></figure>

**Fig. 7a,b.** Effect of a localized velocity change, mimicking a velocity change in a superficial magmatic chamber. Synthetic seismograms are computed with a reflectivity method in the model shown in (a). The velocity is perturbed by 1% in the shaded area and remains constant in the first shallow layer. The CSMW analysis between the perturbed areas and the seismograms is shown in (b). The direct P- and S-waves show no delay: they propagate outside the vertical dyke. The coda exhibits a delay which increases linearly as a function of propagation time reflecting a velocity change in the mimicked magma chamber.
<figure><img src="image_10.png" /><figcaption>**Fig. 8a,b.** Effect of the intrusion of a vertical dyke modeled. Synthetic seismograms are computed with a P-SV finite differences method. a) The model used for the computation: a 30 m thick dyke intrudes in a random medium with 10% velocity fluctuations. The star is the location of the two earthquakes of a doublet and the triangles are the stations for which seismograms are presented in (b). b) Horizontal seismograms recorded at the surface by the four stations shown in (a) and the corresponding CSMW delay diagrams. The black arrows show the theoretical arrival times for a P-to-P diffracted wave at the top of the dyke. The lines in (a) represent the possible depths of P-diffractors for the times of the arrows in stations 1 and 3.</figcaption></figure>

### 3.1.4. A vertical dyke intrusion

In this last case, a vertical dyke intrudes at a depth of 1500 m in the time interval between the occurrence of two earthquakes of a doublet. Four stations located at the surface record the events. This case is simulated using a finite differences computation. A 30 m thick vertical dyke corresponding to a 40% velocity decrease has been inserted in a heterogeneous medium with 10% velocity random fluctuations. Figure 8a,b shows examples of computed seismo-

grams together with their CSMW analysis. Assuming a constant P-wave velocity for waves diffracted by the dyke (no converted waves), one can report the expected time of arrival on the CSMW diagram. The corresponding arrows (fig. 8b) coincide with the occurrence of the first delay discontinuity observed on the diagram. Reciprocally, one could pick up the time of these delay discontinuities at several stations and migrate them. This is done in the model in fig. 8a where the ellipses of possible diffractors for station 1 and 3 are plotted. In our case, the
time uncertainty for picking up the delay dis-
continuity is of the order of the pulse half
width, 0.05 s, giving a 300 m location uncer-
tainty denoted by a dotted zone in fig. 8b.

3.2. *Example of an observation implying temporal changes in velocity inside a volcano*

Figure 9 shows two earthquakes that were recorded before the eruption of Merapi on February 2nd 1992. The first earthquake occurred on 91/06/22 and the second on 91/09/04. The slope of CSMW curve is calculated using a weighted linear regression that takes into account the coherency of the signal in each

window. A decrease in the delay along the
seismogram indicates an increase in velocity as
a function of date before the eruption. The pat-
tern between the direct waves (body waves) is
different from the pattern in the coda: for the
direct waves, within a 2 s interval following
the first arrival, the time delay is constant and
then it gradually decreases along the coda pro-
portionally to the running time. Such data were
used by Ratdomopurbo and Poupinet (1995) to
detect an increase in velocity close to the sum-
mit of Merapi before the February 1992 erup-
tion. This increase in velocity is interpreted as
due to an increase in pressure around the su-
perficial magma chamber starting about 10
months before the eruption.

<figure><img src="image_11.png" /><figcaption>Fig. 9. Example of the CSMW analysis of two seismograms observed before the eruption of the Merapi volcano (Indonesia) on February 2nd 1992. The first events occurred on 91/06/22 and the second on 91/09/04. A S-velocity change of the order of 0.5 per cent beneath the source and the station explains the observed pattern. Such a velocity change may result from an increase in pressure around the superficial magma chamber of Merapi. This increase in pressure would have started about one year before the 1992 eruption of Merapi.</figcaption></figure>
## 4. Conclusions

The analysis of multiplets of micro-earthquakes provides a means to relocate volcanic events precisely and to map deep fractures inside the edifice. For instance, many volcanic events have a low frequency content and it is nearly impossible to locate them with standard techniques because picking up *P*-arrival times in a sufficient number of stations is impossible. We show that these events can be relocated by measuring the difference in arrival times by the cross-spectral technique (precision within 10 m in relative location). Moreover the application of the Cross-Spectral Moving-Window technique to doublets allows *S*-arrival time differences to be picked up very precisely and very small differences to be read in the coda. It can also be used to enhance secondary phases that may be related to the deep structure in the volcano or deep in the crust. We illustrate the use of temporal doublets by computing synthetic seismograms for models that mimic temporal changes of seismic velocity inside the volcano. We analyze these synthetic doublets with the CSMW technique and show that it is possible to measure different types of temporal variations inside a volcano. Several cases are explored:

1) a velocity change in the entire volcano;

2) the appearance of a horizontal dyke;

3) a change in velocity restricted to a circular magma chamber.

We process two events recorded on the Merapi volcano 6 months before the February 1992 eruption. The CSMW technique detects a dependency of $\Delta\tau$ on the arrival time in the coda that is interpreted as a temporal velocity variation in the medium beneath the source and the station. This velocity increase can be related to the building up of pressure inside the magma chamber before the eruption of 1992.

## Acknowledgements

This work was performed within the frame of the French-Indonesian research project on Merapi funded by DRM, Ministry of the Environment and the Ministry of Foreign Affairs, on the French side and by the Volcanological Survey of Indonesia on the Indonesian side.

## REFERENCES

AKI, K. and B. CHOUET (1975): Origin of coda waves: source, attenuation, and scattering effects, *J. Geophys. Res.*, **80**, 3322-3342.

DEICHMANN, N. and M. GARCIA-FERNANDEZ (1992): Rupture geometry from high-precision relative hypocentre locations of microearthquake clusters, *Geophys. J. Int.*, **110**, 501-517.

FREMONT, M.-J. and S. MALONE (1987): High precision relative locations of earthquakes at Mt. St. Helens, *J. Geophys. Res.*, **92**, 10223-10236.

GOT, J.-L., J. FRECHET and F.W. KLEIN (1994): Deep fault plane geometry inferred from multiplet relative relocation beneath the south flank of Kilauea, *J. Geophys. Res.*, **99**, 15375-15386.

LEE, W.H.K. (Editor) (1989): *Toolbox for Seismic Data Acquisition, Processing and Analysis*, IASPEI Software Library, Seismol. Soc. Am., El Cerrito, Ca, vol. 1, pp. 278.

OKADA, H., H. WATANABE, H. YAMASHITA and I. YOKOHAMA (1981): Seismological significance of the 1977-1978 eruptions and the magma intrusion process of Usu volcano, Hokkaido, *J. Volcanol. Geotherm. Res.*, **9**, 311-334.

POUPINET, G., F. GLANGEAUD and P. CÔTE (1982): *P*-time delay measurement of a doublet of microearthquakes, in *Proceedings IEEE ICASSP 82*, Paris, 3-5 May, 1516-1519.

POUPINET, G., W.L. ELLSWORTH and J. FRECHET (1984): Monitoring velocity variations in the crust using earthquake doublets: an application to the Calaveras fault, California, *J. Geophys. Res.*, **89**, 5719-5731.

SCHERBAUM, F. and J. WENDLER (1986): Cross spectral analysis of Swabian Jura (SW Germany) three-component microearthquake recordings, *J. Geophys.*, **60**, 157-166.

VIRIEUX, J. (1986): *P-SV* wave propagation in heterogeneous media: velocity-stress finite-difference method, *Geophysics*, **52**, 889-901.

RATDOMOPURBO, A. and G. POUPINET (1995): Monitoring a temporal change of seismic velocity in a volcano: application to the 1992 eruption of Mt. Merapi (Indonesia), *Geophys. Res. Lett.*, **22**, 775-778.