Introduction
Analysis of seismic data from an array of sensors has been very useful to study wave propagation, particularly, to discriminate between natural earthquakes and underground nuclear explosions. With this purpose, in the 1960s and 1970s, the use of large aperture seismic arrays, generally designed for optimal detection of regional and teleseismic events, was widely spread. An example was, the LASA (Large Aperture Seismic Array) located in Montana (USA), an array of 525 shortperiod seismometers covering a radius of 100 km by groups of 25 sensors concentrically arranged. In each group, the sensors were equally distributed in an area of 3.5 km radius (Green et al., 1965). Another large aperture arrangement was NORSAR, located north of Norway. With a similar configuration, NORSAR consisted of 132 short-period seismometers and covered a radius of 50 km. The sensors of these two arrangements recorded only the vertical component of the movement. Over the years, the use of arrays of seismic sensors of a single component has been very helpful in the first arrival detection and in the identification of apparent velocities and directions of the fields of incident waves. Ringdal and Husebye (1982) made a detailed review of these experiences. Since the first experiences, diverse theories have been proposed aimed to characterize the motion in terms of power density estimators in the ƒ - k domain.
In general, it has been assumed that the motion, recorded by an array, is a stationary random process with zero mean and standard deviation σ. When the motion is coherent, the variance of the process is a good approximation estimator of the power spectral density (Yaglom, 1962). However, the resolution of this conventional approximation estimator is rapidly lost when the coherence of records decreases. As a result, high resolution approximations estimators have received more attention. One of them is the maximum likelihood filter (Capon et al., 1967; Capon, 1969; Kværna and Doornbos, 1986), especially recommended for irregularly distributed arrays.
In the analysis of tele seismic data, it is common to classify the wave field’s components in terms of their apparent velocities (Aki and Richards, 1980). In many cases, the particle motion analysis, as the one proposed by Jurkevics (1988), can be used as a complement in the recognition of wave’s type. Mykkeltveit et al. (1983) note that this procedure allows to distinguish prominent phases in seismograms. The ƒ - k analysis applied to seismic data is also known as Multichannel Analysis of Surface Waves method (MASW; see Park et al., 2007).
The accelerograms of the Michoacán Earthquake of September 19, 1985 (Ms 8.1) recorded by some stations within the Mexico City Valley showed, among many other aspects, the coherence between records (Campillo et al. 1989) and allowed to identify the presence of conspicuous waves of long and short period. This fact served to calibrate a regional crustal model. After this earthquake, the MCAN was deployed. The ensuing strong-motion recorded in Mexico City, clearly show large amplifications beyond the level predicted using simple onedimensional (1D) shear models (Sánchez Sesma et al., 1988; Kawase and Aki, 1989). This strongly suggests significant lateral effects of three-dimensional (3D) nature. The spatial variability and polarization of observed ground motion have been interpreted as 3D effects (Pérez-Rocha et al. 1991; Sánchez-Sesma et al., 1993).
Currently, the MCAN has over 70 three-axial surface accelerometers irregularly distributed in the city; the area covered is of approximately 600 km2. More than 170 intense and moderate events, mainly from the coast, have been recorded by this network.
In this work, seismic records of motion caused by a well recorded earthquake in Mexico City Valley (March 20, 2012), were analysed. The results reveal the propagation of wave components across the Valley in both longitudinal and transversal motions. The characteristics of these components allow identifying the interaction of incoming waves with the geological structure of the Mexico City Valley.
In order to estimate the power spectra of the kinematic wave components in the ƒ - k domain, a wave decomposition was performed, implementing a maximum likelihood filter, based on the Helmholtz theorem and the Capon method (1969). A spectral representation that displays amplitudes, directions and propagation velocities of waves with longitudinal and transversal polarization allows interpreting the behaviour of waves trapped in a sedimentary environment.
Power spectral density of the kinematic components of a wavefield
Fourier analysis in time domain is an usual tool in seismology and it is the keystone of many developments in mathematical-physics. Its extension to space-time configurations in which one or various coordinates are dealt with analytically, shows the great resolving power of Fourier analysis. This is the case when it is used in regular arrays. In many instances it comes out that arrays in which data is avilable, are not regular and even an approximate description is desirable.
An analysis of the kinematic components of a vector wavefield is presented, starting from the continuous description of a vector field in the space-time frame and using the results within the frame of approximate spectral estimation based on the maximum likelihood theory of Capon (1969).
Power spectral density
The displacement field u (x; ω) due to the propagation of a harmonic plane wave is given by:
where u0=u0 ( ω ) is the spectrum of the signal with information of its waveform. Here a = (ax, ay)T = ( ω / c) n is the wavenumber vector with components, ax = ωnx / c and ay = ωny / c; ω is the angular frequency, and c the wave propagation velocity. The wave direction is given by the vector n = (nx, ny)T. Finally, x = (x, y)T is the position vector.
After the spatial Fourier transform, the field in wavenumber domain is:
where δ (⋅) is the Dirac’s Delta function:
Therefore, the spectrum in the wavenumber domain given by eq. (2) is theoretically zero for the coordinates indicated by the vector k = (kx, ky)T that differs from a = (ax, ay)T. This represents a prominent peak in the spectrum and leads to determine the direction and velocity of propagation of a plane wave.
The velocity c and the azimuthal direction of propagation θ are obtained with the following equations:
It is clear that the azimuth is measured with respect to North-South (NS) direction. Capon (1969) designed a high resolution filter to estimate the ƒ- k power spectral density |u(k; ω )|2 when the field is only known in M sensors distributed in the x-y space. This estimation allows the undistorted passage of any monochromatic plane wave that propagates with wavenumber vector a:
where gjl (ω) is the jl element of the matrix G=F−1 and F is the spectral correlation matrix; each element jl of F is given by fjl( ω )=uj (x; ω ) ul*(x; ω ) (here * denotes complex conjugate).
Kinematic decomposition
Assume that P waves (with longitudinal polarization of particle motion) and S waves (with transversal polarization of particle motion) propagate in a medium. According to the kinematic decomposition of Helmholtz, if θ and Ψ are the displacement potentials for P and S waves, respectively, with ∇⋅Ψ = 0, then the total displacements field is written as:
By definition u P=∇Φ and u S=∇×Ψ. Then, it can be written
The total horizontal displacements in the plane x − y, admit the next representation:
since Ψ= (0, 0, Y). If the potentials Φ and Ψ are harmonic and represented by:
Then, the components of the displacements field in the x and y directions acquire the forms:
If l=(ky, − kx)T , the total field u(k; ω ) can be written as:
where k ⋅ l = 0.
In this result, the propagation vectors associated with P and S waves are recognized parallel and perpendicular to the vector k, respectively. Moreover, the displacement field produced by P waves is the component of the total field in the direction of the vector k, whereas the field produced by S waves is the component of the total field in the direction of the vector l:
The unit vectors n and m indicate the directions of propagation of the wavenumbers k and l, respectively, and are given by,
where This kinematic decomposition is developed for the plane x - y, therefore the longitudinal component of Rayleigh waves is analogous to P waves in the above development. In the same way, the transversal component of Love waves is regarded as S waves.
Numerical example on a 2D space
The kinematic decomposition is validated in a 2D setting by simulating the propagation of two plane waves. Thus, two plane P and SV waves are considered with propagation velocities of 2.6 and 1.5 km/s, respectively. The respective azimuths are 60 and -20 degrees. Both waveforms are Ricker wavelets with characteristic frequency of 1 Hz.
The array used in this example is shown in Figure 1. It is composed by a random distribution of 30 stations in a square of 4 km x 4 km. The purpose of applying the ƒ- k spectral analysis in a random array with few stations is to check the strength of the kinematic decomposition formulation to describe the total field in longitudinal and transversal components.
Figure 2 displays the kx - ky diagrams of the longitudinal and transversal displacements, associated with the frequency of 1 Hz. The component distinctions emerge from the application of the separation vector described in eqs. (16) and (17).
In each diagram, energy concentrations can be observed, associated with the incidence or emission of coherent waves observed in a large number of station arrays detected by this analysis. Using the wavenumbers (kx, ky) in which each peak is located, the wave propagation speed c, and direction, θ, can be determined by applying eqs. (5) and (6). The direction corresponds to the azimuth, which is measured from NS direction or in this particular case from axis. Therefore, kx = (c/ ω ) sin θ and ky = (c/ ω ) cos θ
Figure 2 clearly depicts a longitudinal wave located in kx / ω = 0.32 km-1 and ky / ω = 0.2 km-1, which represent a wave traveling with c = 2.63 km/s and θ = 58 degrees. A transversal wave is also observed in kx / ω = -0.22 km-1 and ky / ω = 0.62 km-1, which corresponds to a wave with c = 1.53 km/s and θ = -19 degrees. The description of these two waves from ƒ- k analysis exactly matches the initial parameters of the incident waves.
Numerical example on a 3D alluvial valley
In order to calibrate this kinematic decomposition, the surface field of an elastic 3D alluvial valley on an elastic half-space under the incidence of plane waves (P, SV, SH and Rayleigh) is separated in terms of longitudinal and transversal spectral motion. The reference results were computed with IBEM for the model presented by Sánchez-Sesma and Luzón (1995) for the scattering and diffraction of elastic waves by soft elastic inclusion models of alluvial deposits in an elastic half-space.
This is a closed irregular croissant-type alluvial valley (V) over a half-space (HS), depicted in Figure 3. The valley is 0.25 km depth, and is 2 km of longitude in its longest side. The free surface of both the valley and the half-space are assumed flat. The shear wave speeds are βV = 1 km/s and βHS = 2 km/s; the Poisson ratios are uV = 0.35, uHS = 0.25, and the mass densities satisfy pV = 0.8 pHS. A quality factor of 100 was assumed for both P and S waves inside the valley. The half-space has no internal attenuation.
The array has 20 x 20 stations equally spaced (Δx = 0.0526 km) within a square. It is located inside the croissant valley on free surface as depicted in Figure 3 (b). Note that the right side of the array cross the valley’s edge. The side’s length of this square array is L = 1 km.
Both dimensional parameters: L in the studied direction, and Δx, control the confidence limits of the wavenumbers in the ƒ- k analysis. In an arbitrary array, the minimum wavenumber that should be analysed is: kmin = Δk/2π = 1/L, and the maximum wavenumber, or Nyquist (Nq) wavenumber, that should be studied is kmax = kNq/2π = 1/2 Δx. In this particular case, these two dimensional parameters are: kmin = 1 cycles/km, and kmax = 9.5 cycles/km.
We consider the incidence of a SV wave with θ = -10 degrees, and a Rayleigh wave with θ= 90 degrees, both of them with propagation velocity of 0.5 km/s. For the incoming waves Ricker wavelets with characteristic period equal to 0.33 Hz, was assumed. Synthetic seismograms are depicted in Figure 4.
The ƒ- k analysis in four time windows of 10 s with an overlap of 2 s: 0 - 10 s, 8 - 18 s, 16 - 26 s, and 24 - 34 s (see Figure 4) was applied, in order to follow waves both incident or emitted by the edges of the valley, and identify their type, velocities and directions of propagation.
The windows’ length depends on the desired level of detail in the analysis. Generally, the records’ total duration can be used, but waves with large amplitudes overwhelm coherent waves with low amplitudes in the coda and thus prevent their detection. The window’s resolution depends on the signal frequency, the likely apparent velocity and the array’s spatial dimension. It is recommended that the time window contains at least three wavelengths to allow that the propagating waveform is contained in the selected window, so that their spectrum amplitude and phase are well retrieved.
For each time window we construct kinematic spectra of the amplitudes, velocities and directions of propagation for the longitudinal, transversal and vertical displacements of the wave field. These spectra are constructed from the wavenumbers (kx, ky) of the peaks detected in the kx - ky diagrams for each frequency using eqs. (5) and (6) to compute c and θ, respectively.
The vertical component is used jointly with the longitudinal one to corroborate the possible presence of Rayleigh waves. The propagation velocities and directions of these waves must be equal in both the longitudinal and vertical components, and in the amplitudes, vertical component must be slightly larger, at least for a prevailing fundamental mode.
At the kinematic spectra for the time window 0 - 10 s depicted in Figure 5 (a), using both the longitudinal and vertical components, a prominent Rayleigh wave i can be seen s observed, with frequency between 0.3 and 0.4 Hz, c = 0.4 km/s approximately and θ = 90 degrees. From the transversal component, a S wave of f = 0.35 Hz, c = 0.5 km/s and θ = -20 degrees can be seen. In the same time span, another S wave with ƒ = 0.3 Hz, c = 0.5 km/s and θ = 135 degrees, almost opposite to the previous one, can also be seen .
Figure 5 (b) correspond to time window between 8 and 18 s. It reveals a principal S wave ƒ = 0.35 Hz, c = 0.4 km/s, and θ = -10 degrees. No prominent waves in the longitudinal component can be seen in this time window.
The results depicted in Figure 5 (c) and (d) correspond to time windows 16 - 26 s and 24 - 34 s, respectively, and reveal very interesting patterns of interference of the refracted waves inside the basin and, in some cases, significant emission of waves is observed.
From Figure 5 (c), that covers time window 16 - 26 s, we can see the propagation of initial waves with θ = 0.5 Hz in the three components. In the longitudinal component it travels with c = 0.3 km/s, and θ = 45 degrees. In the transversal component, it revels three S waves with c = 0.3 km/s, in a range of directions of approximately 5, 65 and 130 degrees.
From Figure 5 (d), corresponding to time window 24 - 34 s, we can see in the longitudinal component: the propagation of a wave with a range of frequencies between 0.7 and 0.77 Hz, with c = 0.3 km/s, and q = -145 degrees, i.e., it is traveling backward. With the same frequencies and velocity but with smaller amplitudes there are two waves propagating with q = 17 and 57 degrees, almost opposite to the first one.
There is also a wave with ƒ = 0.4 Hz, with c = 0.35 km/s, and q = 40 degrees. In the transversal component, the propagation of three waves is spotted for an interval of frequencies between 0.7 and 0.77 Hz, with velocities in the interval of 0.27 and 0.32 km/s, and directions of -30, 11 and 125 degrees. Moreover, one can infer the propagation of a transversal cylindrical wave with ƒ = 0.5 Hz, velocity around 0.3 km/s, coming from N135E, i.e., traveling backward.
We attribute the changes in the propagation directions and the presence of several peaks in the power spectra as an evidence of the scattering and diffraction of the incident wave field by the basin’s edge.
Figure 6 shows the kx - ky diagrams associated to only one frequency for each of the four time windows analyzed. The various panels presented give support to the kinematic spectra depicted in Figure 5.
Analysis of seismic data recorded by MCAN
The important site effects that occur at Mexico City subsoil due to earthquakes, especially the long duration and the increased amplitudes on seismograms recorded in free surface inside de City, lead to the search of ways to improve the understanding on the mechanisms that control these phenomena.
In this study, we apply a time-space spectral analysis to obtain in the ƒ- k domain, spectral power densities of the strong ground motion recorded by the MCAN. This method will allow us to figure out the behaviour of the incident and emitted wave fields within the Valley.
Figure 7 shows the distribution of the MCAN as well as the Geotechnical Zonation location. Two parameters control the confidence limits of the wavenumbers in the ƒ- k analysis. They are: the array length in the studied direction, L = 30 km for both sides, and the average separation between stations, Δx = 3.33 km which comes from the maximum instrument density located in the central area of the network. Therefore, taking into account these two values, the minimum wavenumber that should be analysed is kmin = Δk/2π = 0.033 cycles/km, and he maximum wavenumber, or Nyquist (Nq) wavenumber, that should be studied is kmax = 0.15 cycles/km.
Previous studies (Campillo et al., 1989, and Sánchez-Sesma et al., 1993) identified, in the vertical component of the displacements recorded inside Mexico City Valley, crustal surface waves coming from the Mexican subduction zone. Campillo et al. (1989) identified two prominent surface waves with characteristic periods of 0.1 Hz and 0.33 Hz in the records produced by the September 19, 1985, Michoacán earthquake (Mw 8.1) and proposed a cortical model that allowed them to postulate that the long-period waveform of 0.1 Hz corresponds to the fundamental mode of Rayleigh surface wave, and is preceded by a short-period waveform of 0.33 Hz which is associated with the propagation of higher modes of Lg surface waves. Sanchez-Sesma et al. (1993) detected also these two particular waves in the MCAN’s records of the April 25, 1989, earthquake (Mw 6.9) by applying the ƒ- k analysis.
In this study, a more recent event (with characteristics similar to the ones mentioned above) was selected to be analysed, i.e., large magnitude with its epicenter in the Mexican subduction zone, and with important site effects when arriving to Mexico City Valley. The event is the March 20, 2012, earthquake (Mw 7.4) recorded at 60 MCAN stations.
The displacements in six MCAN arbitrary sites due to this event after the double integration of the accelerograms, appear in Figure 8. In the vertical displacement of the set of records, it is also possible to observe two conspicuous incidences: Rayleigh (ƒ = 0.1 Hz) and Lg (ƒ = 0.33 Hz) surface waves between 80 and 130 s.
A successive ƒ- k analysis of time windows can determine the wave field behaviour of the records at late times, which are very likely the result of multi-pathing inside the Valley, i.e., the product of waves interaction with the edges of the basin. This interpretation should be subjected to closer scrutiny considering more realistic cortical models and the correlation with surficial geological structures.
In this sense, the displacement spectral densities and the kinematic decomposition in longitudinal and transversal components, for four time windows: (a) 80 - 130 s, (b) 120 - 170 s, (c) 160 - 210 s, and (d) 200 - 250 s were estimated. The last three time windows correspond to the early coda of the earthquakes.
Figure 9 displays the kinematic velocity spectra, direction of propagation and Fourier amplitudes for the four time windows in the frequency band between 0.06 and 0.6 Hz. These spectra were constructed with the three maximum peaks of the diagrams for each frequency. This representation allows identifying the spectral characteristics of the kx - ky wave field components.
Figure 9 (a) shows three representative peaks in all components at frequencies equal to 0.11, 0.15 and 0.33 Hz. The first peak with relevant amplitudes can be observed in an interval of frequencies between 0.1 and 0.13 Hz which propagates in direction around -10.3 degrees, with velocity equal to 2.9 km/s. This peak appears in the longitudinal and vertical direction, which strongly suggests that it is the fundamental mode of Rayleigh waves. Another representative peak occurs at 0.33 Hz in the transversal component, which propagates at -13 degrees with a velocity of 3.6 km/s. We associate this wave with Lg surface waves. The direction of propagation of these two waves is very close to the epicentral direction N10W. This fact suggests that they are incident waves coming from the source. The inferred type of wave and velocities matches with the ones identified by Campillo et al. (1989). The speed and direction of propagation of these two wave episodes are directly seen on seismograms. Therefore, we consider it as a validation of the method that can track waves crossing the Valley and are not clearly detected on seismograms. Another identified waves episode in the transversal component corresponds to the ƒ = 0.15 Hz, with c = 2.8 km/s, and θ = -18 degrees.
Figure 9 (b), (c) and (d) which correspond to the early coda of the seismograms show the kinematic spectra of the longitudinal, transversal, and vertical components. It can be seen, in the interval of frequencies between 0.2 to 0.5 Hz, comparable amplitudes in the three kinematic components with significant changes in directions and velocities of propagation. As in the synthetic model case of a croissant-type alluvial valley, we attribute this behavior of the wave field to an evidence of the diffraction of the incident field due to its interaction with the edges of the basin, and in this real case, also to the multipath that waves take from the epicenter to the Mexico Valley. Moreover, for these frequencies our results suggest that the field is composed by surface waves higher modes, Rg and Lg, in agreement with Campillo et al. (1989).
The phase velocities obtained from the analysis of the MCAN data, in the range of frequencies between 0.1 to 0.5 Hz are consistent with those obtained from the dispersion curves of the accepted theoretical crustal models for central Mexico (see Gaite et al. 2012). A discontinuity in the experimental results about ƒ = 0.2 Hz, both in phase velocities and propagation directions, suggests a possible conversion of surface wave to higher Rg and Lg modes, probably generated on the south edge of the Transmexican Volcanic Belt. Apparently, the cortical paths from this geological structure towards the Valley of Mexico, is very efficient for the propagation of these mode shapes. This follows from the conspicuous amplitudes of the power spectrum of these forms of propagation and could be related with the regional amplification identified by Ordaz and Singh (1992).
The resolution of the power density estimators decreases rapidly for frequencies larger than about 0.5 Hz. This is a consequence of the decrease of the signal-to-noise ratio mainly because the minimum distance between MCAN sensors is around 3 km and wavelengths are less than 6 km, because waves propagate typically with horizontal velocities less than 3 km/s.
Conclusions
The importance of estimating the ƒ- k power spectral densities within the analysis of wave propagation from data of an array of sensors has been discussed. Higher resolution filters based on Capon’s (1969) maximum likelihood estimators, were applied. We have also proposed a kinematic decomposition of the wave field in longitudinal and transversal components, which permits to observe the directions, velocities, and polarizations of the surface displacement field.
For validation, the kinematical decomposition was implemented for two synthetic models, in 2D and 3D configurations, respectively. The results were successful because the parameters of the incident wave field like wave type, velocity and direction of propagation were well identified. The 3D model was a croissant-type alluvial valley, in which the influence of basin geometry was identified in surface motion.
For the March 20, 2012, earthquake recorded by the MCAN, the analysis and the kinematic decomposition in longitudinal and transversal components permitted to detect the propagation of waves in different time windows. In early times of seismograms, the incident wave field is recognized, while for late times, the analysis made possible to follow the variety of waves with frequencies between 0.2 and 0.5 Hz that are propagating across the Valley. They may be the result of the interactions of the incident wave field with the edges of the basin. The analysis of early times also reveals locally generated surface waves. They are somehow shadowed by incident waves, but their distinct azimuths may allow their detection. These facts may clarify the retrieved surface waves from correlations of strong ground motion for pairs of observers inside the Valley.
Acknowledgements
Thanks are given to A. Arciniega-Ceballos for her suggestions and detailed review of the manuscript, and to M. A. Contreras-Zazueta for computing the frequency response for the array within the croissant-type alluvial valley. The help of E. Plata and G. Sánchez and their team of Unidad de Servicios de Información (USI) of Instituto de Ingeniería UNAM was crucial to locate useful references. This work was partially supported by AXA Research Fund, and by Dirección General de Asuntos del Personal Académico (DGAPA) of UNAM, under Grant Number IN104712.