## Abstract

The zero-lag cross-correlation technique, used for array analysis in the hypothesis of plane waves, has been modified to allow the wave front to be circular. Synthetic tests have been performed to check the capability of the method, which returns the input test data when the source–array distances are not greater than two or three times the array aperture. For this distance range the method furnishes an estimate of the apparent velocity and the epicentral coordinates of the source. For more distant sources the method becomes equivalent to that based on the planar-wave approximation, which gives an estimate of the backazimuth to the source and the apparent velocity. The method has been applied to seismic data recorded at the active volcano located at Deception Island, Antarctica. 35 volcanic long-period events occurring in a small swarm were selected. Results show that the epicentres are close to the array (between 0.4 and 2 km) and aligned in a SW direction, in agreement with one of the main directions of the fracture system of Deception volcano.

## Introduction

An array, or ‘seismic antenna’ (Chouet 1996a), is a network of seismic stations which can record coherent wave-packets. The characteristics of the arrays are determined by the nature of the seismic waves and the distance to the source. One of the main applications of large-aperture seismic arrays has been the discrimination and control of nuclear explosions (e.g. Mykkeltvert *et al.* 1990), but they have been extended to other applications such as the analysis of regional activity (Vogfjord & Langston 1990), the analysis of seismic scattering (Dainty & Toksöz 1990; Del Pezzo *et al.* 1997), the study of the scaling laws of body-wave spectra (Berger *et al.* 1984), source parameters (Fletcher *et al.* 1987), measurements of surface-wave dispersion curves (Horike 1985) and the monitoring of the seismic activity in volcanic areas (Neuberg *et al.* 1994). In volcanic seismology, small-aperture, short-period seismic arrays have been used in different areas, e.g. Kilauea (Goldstein & Chouet 1994), Stromboli (Saccorotti *et al.* 1998), Vesuvio (De Luca *et al.* 1997), Teide (Del Pezzo, Ibáñez & La Rocca 1997) and Deception Island (Almendros *et al.* 1997), to monitor the seismic activity. This last application differs from the above basically in two aspects: (1) the distance between the source and the instruments is short, and (2) the nature of some of the signals is different from tectonic earthquakes. The seismic arrays are often used to track the seismic source in cases of multiple or complex volcanic sources, in order to focus on the source dynamics. Other applications have been the characterization of surface wave types (Ferrazzini, Aki & Chouet 1991) and the location of the volcanic tremor source (Chouet *et al.* 1997). This last one is important in volcanic seismology because it could be the only way to locate quakes such as long-period events or volcanic tremor, which cannot be easily located by picking arrival times.

Array techniques can be developed in the frequency or time domains to locate coherent wave-packets. All frequency-domain methods are based on the measure of the cross-spectral matrix, which contains the cross-spectrum evaluated at all the pairs of stations that constitute the array. For example, the MUSIC algorithm (Goldstein & Archuleta 1991) has been applied to study the wave field of volcanic quakes at Kilauea (Goldstein & Chouet 1994) and Stromboli (Chouet *et al.* 1997). Time-domain methods are based on the calculation of the array-averaged cross-correlation between all the station pairs as a measure of the coherence of the wave-packet across the array. Among them, the zero-lag cross-correlation method (Frankel *et al.* 1991) has been applied in other volcanic areas such as Teide volcano (Del Pezzo *et al.* 1997) and Deception Island (Almendros *et al.* 1997). The array techniques (both in frequency and time domains) are based on the plane wave approximation, which is realistic for epicentral distances greater than approximately four or five times the array aperture. For closer distances it becomes necessary to modify this model taking into account the wave-front curvature. In seismological monitoring of active volcanoes, owing to the low magnitude of the volcanic quakes and the high attenuation of the volcanic rocks, it is often necessary to deploy the instruments as close to the volcanic quake source as possible, sometimes at distances shorter than 1 km (e.g. Goldstein & Chouet 1994, at Kilauea Volcano).

In the present work we introduce the circular-wave-front geometry into the zero-lag cross-correlation array method. We have chosen this time-domain method because: (1) it permits us to track the evolution of the volcanic source when the signal duration is short or when the delay between two coherent wave-packet arrivals from the same source is small; (2) the results are less sensitive to the window duration than in the case of the frequency-domain methods; and (3) the circular-wave-front geometry can be easily implemented in the method. We will show that, for short epicentral distances, the plane-wave-front model leads to biased solutions, while the circular-wave-front technique gives more exact results. Moreover, the use of the circular-wave-front hypothesis allows us to estimate the epicentre location. This determination is important when a clear first onset is not available (for example, volcanic quakes such as tremor or long-period events). We applied the circular-wave-front technique both to synthetic data and to real seismograms belonging to a subset of volcanic quakes recorded in an experiment carried out at Deception Island (Antarctica). We show that it is possible to estimate the epicentral position and the apparent slowness for nearby sources, obtaining more information than the plane wave technique can provide.

## Method

The zero-lag cross-correlation method (Frankel *et al.* 1991) has been used to analyse several types of seismic signals (see e.g. Mori *et al.* 1994; Del Pezzo *et al.* 1997; Almendros *et al.* 1997). It is based on the evaluation of the cross-correlation between signal pairs. It can be applied to very-short-duration signals, because the results are not as strictly dependent on the time window length as are the methods based on FFT. In fact, a window that contains just one cycle of the signal is enough to locate the wave packet. The main assumptions that have to be taken into account are: (1) there are not different local effects at the sites of the array sensors; (2) the seismic stations are located in the same plane (ideally horizontal); and (3) the wave fronts that propagate across the array are plane. These conditions are satisfied if the array site is geologically homogeneous, the topography is smooth and the seismic source is far enough from the array. In this case the zero-lag cross-correlation method supplies information about the apparent velocity of the waves and the backazimuth to the source. To avoid the restriction concerning the distance between the source and the array, and therefore, to be able to apply this method to nearby sources, the plane-wave-front hypothesis needs to be replaced.

### Formulation of the problem

Let the array be composed of *N* sensors at locations **r**_{j} ( *j* = 1 … *N*) with respect to a reference point (Fig. 1). The sensors sample the wavefield in both time and space, providing the seismic traces *W*_{j} (*t*) ( *j* = 1 … *N*).

Let us suppose that the spatial and temporal properties of the wave front generated by a unique source can be characterized by a set of *m* independent parameters {**α**_{i} , *i* = 1 … *m*}. Our problem is to estimate those parameters, which in terms of matrices can be expressed as the one-column matrix **α**.

The cross-correlation function is defined by

where *W*_{j} (*t*) is the seismic trace at the station *j* and 〈〉 indicates the time average. Let *τ*_{jk} be the time delay between stations *j* and *k*. *R*_{jk} will reach a maximum when the traces are shifted by *τ*_{jk} . The normalized average of *R*_{jk} over all the possible pairs is

This function *R* will take the maximum when all the trace couples are shifted by *τ*_{jk} . If *τ*_{jk} can be expressed as a function of **α**, the problem of finding the parameters that return the maximum averaged zero-lag cross-correlation can be solved.

#### Plane waves

Let us suppose that the seismic source is located in a homogeneous half-space, so far from the array site that we can consider the incoming wave fronts to be planar (Fig. 1a). In these conditions, the source is at infinity and the phase delay between the array stations is dependent only on the apparent slowness vector, **S** ≡ (*S*_{x} , *S*_{y} ). Alternatively, in polar coordinates we can use the backazimuth to the source, *A*, clockwise from north, and the apparent slowness-vector magnitude, *S*, related to the previous one by the expressions

The expected arrival time to the *k*th station can be expressed as

where *t*_{0} is the arrival time to the origin of the coordinate system. Then, the delay between stations *j* and *k* is given by

Therefore, in the case of plane waves the matrix **α** becomes (*S*_{x} , *S*_{y} ) in cartesian or, alternatively, (*S*, *A*) in polar coordinates:

#### pherical waves

If the distance to the source is not large enough compared to the array size, the assumption of plane wave fronts is not justified, and it is necessary to consider spherical wave fronts (Fig. 1b). Hypothesizing a seismic source located at coordinates (*x*_{s} , *y*_{s} , *z*_{s} ), the waves generated at the source propagate with radial apparent slowness *S* outward from the epicentre. If we use cylindrical coordinates, with the origin at the epicentre, this apparent slowness field can be expressed as

where *S*_{0} is the actual slowness for the waves in the medium and *ρ* is the epicentral distance. The time delay between stations *j* and *k* is given by

For spherical waves in a homogeneous medium the matrix **α** has four elements: the three source coordinates and the actual slowness. In the case of a shallow source the number of parameters reduces to three. Therefore, if we assume a point-source located at the surface the components of matrix **α** are the two epicentral coordinates and the slowness magnitude. In polar coordinates, with the previous notations, we can use as parameters *S*, *A* and the epicentral distance *D*:

The waves spend a time *t*_{k} to reach the *k*th station starting from the source, given by

and the time delay between stations is

For a source at depth within a vertically inhomogeneous medium, the intersection of the wave front with the surface will be circular, and therefore expression (12) remains valid.

### Solution of the inverse problem

The inverse problem can be solved in general by using a grid search method in the parameter space. This search works as follows: each parameter *α*_{i} is ranged over an interval [*α*^{mini} , *α*^{maxi }] in which we are sure that the solution lies. The grid spacing, Δ*α*_{i} , has to be chosen to be small enough to improve the resolution of the method, but a limitation comes from the required computer time. The number of steps for each parameter is

Therefore, we are working over a *q*_{1}* × q*_{2}* × … × q*_{m}*m*-dimensional grid. Each point in the grid is defined by the indices (*n*_{1} , *n*_{2} , … , *n*_{m} ), in such a way that the corresponding parameters are

Several sets of predicted arrival-time delays between the *N* stations of the array, {*τ*_{jk}, *j = 1 … N*, *k = j + 1 … N*}, are generated for each point in the grid, using, for example, eqs (1) or (2). Then, the seismograms are shifted to align them in time and the array-averaged zero-lag cross-correlation is calculated. Its maximum, which we call maximum averaged cross-correlation (hereafter, MACC) will correspond to the best estimate of *α*, and therefore to the solution.

In practice we have to ensure that the grid amplitude is wide enough to include the solution. For example, in the case of polar coordinates, we have to ensure that the grid search for the parameter *A* occurs between 0° and 360°. The useful range for the slowness depends on the kind of waves we want to study. In the case of surface waves in volcanic rocks, it may lie in the 0.5–2.5 s km^{−1} range. The time delays between stations and the array-averaged cross-correlations need to be evaluated for all the grid points. For a 3-D grid (circular-wave-front geometry), it is computationally quite hard. We applied first the plane wave approximation to obtain a crude estimate of the backazimuth and the slowness, and then used the 3-D grid search, letting the epicentral distance be fully variable, and constraining the backazimuth and slowness to vary only around the first solution.

## Synthetic Tests

Hereafter we refer to CWM (circular-wave-front method) and PWM (plane-wave-front method) as the application of the zero-lag cross-correlation method using the circular and the plane-wave-front approximation, respectively.

To study the capabilities of CWM a synthetic test has been performed. The synthetic signals are supposed to be originated by an isotropic shallow source, and recorded at the 18 seismic stations of a small-aperture seismic array that simulates the one deployed in Deception Island, South Shetland Islands, Antarctica (Fig. 2) during the 1994–1997 austral summer surveys (Ibáñez *et al.* 1997). The method should ensure that the input parameters are well reproduced. The test was developed to compare CWM with PWM and to check the limitations imposed by the size of the array and the effect of local irregularities of the wave front. The shape of the test signal is designed to simulate the long-period events recorded at Deception Island, whose main characteristics are (Del Pezzo *et al.* 1998): (1) quasi-monochromatic spectral content, with dominant frequencies around 2 Hz; (2) lack of impulsive phase arrivals, which makes the location difficult with techniques different from those associated with arrays; (3) spindle-shaped envelope; and (4) they propagate with apparent slowness around 1.4–1.5 s/km (they are composed mainly of surface waves).

We used a signal shown in Fig. 2(b) and (c), given by

The values of the constants are *A* = 100 m s^{−1}, *B* = 4, *t*_{0} = 0.1 s, * f*_{0} = 2 Hz. The case *B* = 1 is just an Ohnaka pulse, but we use a higher *B* value to obtain a less impulsive arrival. This pulse is supposed to propagate with circular wave fronts and apparent slowness of 1.4 s km^{−1} from 378 single points on the surface, distributed over a variety of backazimuths and distances (Fig. 2a) given by

Backazimuths cover the whole horizontal plane and distances range from 0.1 to 9 km. The slowness is fixed at a value consistent with the propagation of surface waves in a low-velocity medium. As we discuss further, it is a realistic approximation in many cases. The signals are contaminated with real noise (signal-to-noise ratio of 10) obtained from noise samples at the array stations. The geometrical spreading has been taken into account by reducing the signal amplitudes by a factor of 1/√*r*. Fig. 2 shows an example of the synthetic traces at the array stations. The parameters for the analysis have been the following: window length of 1 s from the beginning of the signal; slowness search grid centred in the actual solution with 0.04 s km^{−1} spacing and a width of 3.2 s km^{−1}; distance spacing for CWM of 25 m, from the array centre to 10 km.

### Source location

We used PWM as well as CWM to estimate the slowness, backazimuth and distance (CWM) for the synthetic events. For each one we applied first PWM and found backazimuth and slowness. Then, we applied CWM in a grid centred in this last solution. In Fig. 3 we report the obtained solution for both PWM and CWM in two distance ranges: 1–9 km (far sources, Fig. 3a) and 0.1–1 km (near sources, Fig. 3b). The error for each solution is reported using different grey scales.

For distant sources, both methods give the same solution for backazimuth and slowness, as could be expected. For near sources, CWM returns the true backazimuth and slowness with an error less than 3° and 5 per cent, respectively. The percentage error has been determined subtracting the estimate from the true value and taking the absolute value. Moreover, this method provides a good estimate of the true distance inside of a percentage error smaller than 20 per cent, up to a distance of 1.5 km from the source. PWM returns biased estimates of backazimuth and slowness with an uncertainty greater than 10° in backazimuth, and greater than 15 per cent in slowness, for sources closer than around 400 m. From Fig. 3b, it can be observed that there is no symmetry in the reported distribution of the difference between the estimated and the true value. This asymmetry is due to the array configuration that conditions the resolution of the methods in backazimuth and slowness.

### Array size

We tested the array size effects using CWM for three different arrays with the same geometry but with different sizes: half, original and double. The aim of this test is to study the possible relation between the distance estimate capabilities and the array sizes. For this experiment, we have fixed the backazimuth (200°N) and the slowness (1.4 s km^{−1} ), and we have ranged the distance of the source from 0.1 to 9 km with the same distance spacing shown in eq. (16). All three arrays locate the source in backazimuth and estimate the slowness of the waves without bias. However, the differences among them are in the uncertainty interval, which is greater for the smaller array than for the largest. It is noteworthy that all sizes provide good results in backazimuth and slowness because the apparent slowness selected was big enough to guarantee that the small array could correctly sample the wavefield. The maximum distance that can be estimated is related to the array size. In Fig. 4 we plot an example of the MACC values versus distance for two synthetic sources located at 0.5 and 1.8 km from the array centre. For the nearby source (solid lines), a clear maximum appears for the two larger arrays, while for the smaller one the maximum can be roughly observed. For the far source (dashed lines), only the largest array can give a coarse solution. We have observed that the critical distance, from which we cannot distinguish a maximum, is around two times the array aperture. This limit is directly related to the time that the incoming wave front spends in crossing the array, and therefore, to the apparent slowness of the waves: the greater the slowness, the better-constrained the distance. From the above results, it could be assumed that the easiest way to improve this method is to design an array as large as possible. But the array aperture is limited by the lack of coherency of the seismic signals as they propagate at long distances compared with the wavelengths involved.

### Wave-front deformations

In order to take into account possible deformations of the circular wave front due to the presence of small lateral heterogeneities between the seismic source and the array site, we have introduced a random error, uniformly distributed, in the synthetic time delays. This random error ranges from 1 per cent up to 50 per cent of the traveltime. The results reported in Table 1 show that, up to an error of 30 per cent, the CWM returns unbiased estimates for a near source located 400 m away from the array centre. Moreover the MACC decreases as the introduced random error increases. Lower values of MACC imply less accurate solutions.

## Data Analysis and Results

Almendros *et al.* (1997) and Del Pezzo *et al.* (1998) have studied the volcanic quakes recorded in Deception Island (Antarctica) at the dense seismic array reported in Fig. 2. They identified, on the basis of their time-domain appearance and spectral content, at least three groups of events: volcanic tremor, hybrids and long-period events. We do not take into account tremor and hybrids, because they are composed of a complex mixture of body and surface waves (Almendros *et al.* 1997). We consider only the long-period events that show an almost pure surface-wave composition, from the first onset to the end of the coda. On the basis of the best signal-to-noise ratio, we selected 35 long-period events belonging to a 42 hr swarm which occurred at Deception Island from February 5 to 7, 1996. In Fig. 5 we show a three-hour recording of part of this swarm. Small variations of the spectral content of the quakes were observed, with the maximum spectral peak varying from 1.4 Hz to 2.5 Hz. The maximum occurrence rate was about 100 quakes per hour. Del Pezzo *et al.* (1998) applied PWM to this type of event using a short time window sliding along the seismogram. This analysis shows that the apparent velocity and the backazimuth are almost constant during the whole signals. We selected for the present analysis a time window of 1 s. PWM was applied to the whole seismogram. Then, CWM was applied to the window that shows the highest MACC value.

The grid search was performed from −3.2 to 3.2 s km^{−1} with a spacing of 0.04 s km^{−1} in *S*_{x} and *S*_{y} . Once the PWM solutions were obtained a new grid search was started for the CWM analysis. The new grid was centred at the best PWM solution with a search interval of ±1.6 s km^{−1} with a spacing of 0.04 s km^{−1} in *S*_{x} and *S*_{y} . The distance spacing was 25 m from the array centre to a maximum distance of 4 km. In Fig. 6 we show an example of a CWM solution. The maximum peak of the array-averaged zero-lag cross-correlation is clearly visible in this case at a distance of 0.4–0.5 km from the array centre. This peak is 20 per cent higher than the MACC value at large distances. There are two main sources of error: one due to the array geometry, time sampling and grid search, and the other due to the seismic noise and the lack of coherence of the signal across the array caused by the propagation factors (ray paths, reflections, etc.). For a more complete description of the error evaluation see Del Pezzo *et al.* (1997). The uncertainty value on the MACC peak (Δ*C*) defines an uncertainty region in the parameter space (Fig. 7) from which errors in backazimuth, slowness and distance can be deduced. The distance error interval is not symmetric, owing to the distribution of MACC that becomes flatter at far distances.

The solutions obtained for the whole data set are reported in Table 2. The percentage MACC improvement of the results obtained using CWM with respect to PWM is also shown. It is noteworthy that backazimuths and slowness values are similar for both methods. Moreover, there are many events for which the improvement is greater than 15 per cent, reaching a maximum of 28 per cent. Twelve solutions show no clear improvement (<5 per cent).

There is a direct relationship between the MACC solutions and the predominant frequency of the quakes. In Fig. 8 we plot two extreme examples. The long-period event that has a MACC improvement of 28 per cent is located at short epicentral distance (400 m); the circular-wave-front approximation is much better than the plane wave-front geometry. In the example with no improvement (1 per cent) we are in the limit of resolution, and the error area is larger than the previous one. On the other hand, the spectral content is different: the closer event shows a frequency peak around 2.8 Hz, while the further one is peaked at a lower frequency, 1.5 Hz. This relationship between the distance and the frequency peak is an indirect check of the goodness of our distance estimate: highest-frequency peak means lower attenuation, due to the shorter path. Such a clear attenuation effect is due to the low-*Q* materials composing the volcanic rocks at Deception Island (Vila, Correig & Martí 1995). In Fig. 9 the obtained solutions are reported together with the error area. Epicentres delineate a source region clearly aligned along 210°N.

## Discussion and Conclusions

We present an improvement of the array method to locate coherent seismic signals based on the cross-correlation analysis in the time domain, introducing the circular geometry for the wave front incoming to the array. We have found that using this improvement (CWM), we are able to locate the epicentre of the seismic source for epicentral distances smaller than about three times the array aperture. This method is useful in volcano seismology analyses, because of the following facts: (1) in many cases we have to deal with the analysis of coherent signals with no clear onset or phases; and (2) owing to the strong attenuation of the volcanic material, the instruments have to be located very close to the source in order to preserve the information on the spectral shape of the source signal.

The method was tested with a numerical simulation of a seismic signal produced by an isotropic point source located at the surface. An explosive source generated by a pressure step occurring in a fluid medium can be regarded as an isotropic point-source, similar to those that we can expect at Deception Island. This type of source is often invoked to explain the occurrence of the so-called explosion-quakes that occur in active volcanoes (e.g. Crosson & Bame 1985). We generated a wavelet of 1 s duration with an emergent first arrival characterized by slowly increasing amplitude. The results of this synthetic test indicate that CWM works better than PWM for short epicentral distances (less than three times the array aperture). PWM gives biased solutions for short epicentral distances, this bias depending on the geometry of the array. However, the most important advantage of CWM is the fact that it provides information on the epicentral distance of the source.

We applied the method also to real data. The selected data come from a set of thousands of volcanic earthquakes recorded at the Deception Island seismic array. We selected a swarm of low-frequency quakes composed predominantly of surface waves. We have obtained solutions for the 35 analysed long-period events. These solutions are not only the estimates of backazimuth and apparent slowness, as PWM provides, but also the estimates of the epicentral distance. In all cases the MACC values are improved by using CWM, in comparison to PWM. The level of improvement is directly related to the distance to the source: the shorter the distance, the greater the MACC improvement. The obtained distances range from 0.3 to 2.5 km, with a clustering of the activity around 0.5 km. This result is confirmed by the spectral content of the signals. The events closer to the array have the spectral peak at a higher frequency than the more distant ones. We interpreted this spectral shift as due to an attenuation effect along a longer path. The strong improvement of the MACC value for short epicentral distances indicates that the hypothesis of circular geometry for the incoming wave fronts is more realistic than the plane geometry.

At short epicentral distance, near-field effects could affect the wave-front geometry. An estimate of the wavelength can be obtained by supposing an apparent velocity of 0.5 km s^{−1} for the waves at 2 Hz, which gives a result of around 250 m. If we assume that the source of the events is modelled by a pure pressure step inside a fluid in spherical symmetry (e.g. Crosson & Bame 1985), the size of the source would be less than 100 m. The observed wavelength is greater than twice the expected source size and little near-source effect would be expected (Madariaga 1989). Therefore, the circular-wave-front geometry could be a good approximation of the real situation, even though we can modify very easily the proposed method to include any other front shape.

The surface-wave composition of the quakes used in the present analysis can be explained with an explosive shallow source, producing body waves (which attenuate strongly and arrive at the array with amplitudes at the same level as the noise) and high-amplitude surface waves. A similar wave composition has been observed for the explosion quakes recorded at the explosive, shallow volcanic source of Stromboli volcano (Neuberg *et al.* 1994; Chouet *et al.* 1997) by using two different arrays located close to the craters. These shocks, except for the first two seconds from the first onset, show a predominant composition of surface waves. It is noteworthy that the epicentres of the events analysed in the present paper are aligned with 210°N. This direction is parallel to a fracture system reported by Martí & Baraldo (1990), responsible for part of the last eruptions that took place on the island and the present fumarolic activity. We hypothesized that a possible source model would be some kind of interaction between water and hot materials, producing a sudden phase change from water to vapour. This would happen at the contact between the shallow aquifers which permeate the whole island (Martini & Giannini 1988) and the underlying hot zone. This model, proposed to explain the source of the volcanic events at Deception Island (Del Pezzo *et al.* 1998), does not correspond to the widely used model for the long-period events (Chouet 1996b), though the final waveform is actually very similar.

As a summary, we have found that for the analysed events: (1) the array delays are better fitted by a circular-wave-front geometry than a planar geometry; (2) the apparent slowness of the waves is high, and therefore consistent with the arrival of surface-wave packets; (3) the spectral content is related to the distance; (4) the epicentres are aligned with a well-known fracture system; and (5) the obtained source region is contained on the island itself. Moreover, our source model is based on these facts: (1) there are no active volcanic vents in Deception Island; (2) the source region is located in an active region, between two areas of historical eruptions (1842 and 1919); (3) there is a shallow aquifer in the zone that, during the survey period, provided an important contribution of thaw water from nearby glacial areas.

Finally, the proposed method can be easily applied to other volcanic signals such as tremor; the only requirement is that the array is close enough to the source.

### Acknowledgements

This paper has been partially supported by the project ANT95-0994-C03-02 and by the Grupo de Investigación en Geofísica de la J.A., 4057. We gratefully acknowledge the Spanish Army and Navy for providing us with logistical help during the Antarctic surveys in the Gabriel de Castilla Station.

## References

*S*wave velocity structure down the basement in urbanized areas

*P*and

*S*wave-fronts from the dense three-dimensional array at Garni, Armenia