gji%2Fggx426

код для вставки на сайт или в блог

ссылки на документ

Confidential manuscript to be submitted to JGI
Site Characterization at Groningen Gas Field Area Through
Joint Surface-Borehole H/V Analysis
Zack J. Spica1 , Mathieu Perton2 , Nori Nakata3 , Xin Liu1 , Gregory C. Beroza1 .
1 Department of Geophysics, 397 Panama Mall, Stanford University, Stanford, California 94305-2215, USA.
2 CONACYT, Instituto de Geofísica UNAM, Unidad Michoacán, Antigua Carretera a Pátzcuaro 8701, 58190, Morelia,
3 ConocoPhillips
Michoacán, Mexico.
School of Geology and Geophysics, University of Oklahoma, Norman, USA.
Key Points:
•
•
•
•
•
Site characterization
Interferometry
Ambient Seismic Field
H/V spectral ratio
Groningen
Corresponding author: Zack Spica, zspica@stanford.edu
–1–
Confidential manuscript to be submitted to JGI
Abstract
A new interpretation of the H/V spectral ratio in terms of the Diffuse Field Assumption
(DFA) has fueled a resurgence of interest in that approach. The DFA links H/V measurements to Green’s function retrieval through auto-correlation of the ambient seismic field.
This naturally allows for estimation of layered velocity structure. In this contribution, we
further explore the potential of H/V analysis. Our study is facilitated by a distributed array
of surface and co-located borehole stations deployed at multiple depths, and by detailed
prior information on velocity structure that is available due to development of the Groningen gas field. We use the vertical distribution of H/V spectra recorded at discrete depths
inside boreholes to obtain shear wave velocity models of the shallow subsurface. We combine both joint H/V inversion and borehole interferometry to reduce the non-uniqueness of
the problem and to allow faster convergence toward a reliable velocity model. The good
agreement between our results and velocity models from an independent study validates
the methodology, demonstrates the power of the method, but more importantly provides
further constraints on the shallow velocity structure, which is an essential component of
integrated hazard assessment in the area.
1 Introduction
Due to long-term exploitation of a large onshore gas field and subsequent compaction of the reservoir at depth, the Groningen area in the northern Netherlands is subject
to ground subsidence and induced earthquakes. To date, the largest induced earthquake
occurred near Huizinge in August 2012 and was recorded with a local magnitude ML of
3.6 (moment magnitude M = 3.4). This event revived concerns about the hazard from induced seismicity and prompted the government and the Nederlandse Aardolie Maatschappij (NAM) to take action. As part of the response, NAM built a dense permanent borehole
microseismic network covering an area of about 35 × 45 km2 . The network (Fig. 1) is
composed of ∼70 accelerographs at the surface that are co-located with a 200-meter depth
borehole in which geophones are installed at depth intervals of 50 m. The network is fully
operational and has recorded continuous waveforms since late 2015.
It is well known that near-surface lithology can strongly influence damage from
earthquake shaking by increasing the amplitude and duration of shaking, and by responding nonlinearly to incident seismic waves [e.g., Olsen, 2000]. For this reason, site characterization is of great importance to seismic hazard analysis [Bommer et al., 2017]. Such
characterization requires a knowledge of geomechanical properties of the stratigraphy at a
site, which can be often challenging to obtain.
Over the last few decades, ambient seismic field seismology has emerged as a valuable tool to characterize shear wave (VS ) velocity models over all distance scales. For site
characterization, these techniques include the inversion of dispersion curves obtained by
ambient seismic field correlation using small-aperture arrays such as the spatial autocorrelation [SPAC; Aki, 1957] or frequency-wavenumber methods [f-k; Capon, 1969; Lacoss
et al., 1969]. Depending on the aperture of the array, these techniques can recover the local 1-D velocity structure from surface to a depth of up to several hundred meters. Ambient seismic field cross-correlation tomography with high-density arrays is an emerging
method to obtain high-resolution 3D models of the shallow surface [e.g., Nakata et al.,
2015; Roux et al., 2016]; however, these kinds of experiments are very expensive, and require significant investment for design, deployment, and operation.
The ratio of the horizontal to vertical (H/V) components of the ambient seismic field
(sometimes referred to as microtremor) is a widely used method to determine a simple
velocity model of the subsurface (i.e., one layer over a half-space) and from that the frequency dependent site response. The reason this method gives reliable results, however,
has long been controversial. This was due to the absence of a clear theoretical basis for
the measurements, which leads, inevitably, to lack of clarity in its interpretation. It is not
even clear which waves comprise the noise field that generates the H/V peak frequencies
[Nakamura, 2000; Malischewsky and Scherbaum, 2004; Bonnefoy-Claudet et al., 2008].
–2–
Confidential manuscript to be submitted to JGI
The possibility to retrieve the elastodynamic Green’s function between two stations
embedded in an elastic medium from the average time-domain cross-correlation of ambient field records [Weaver, 2005] led Sánchez-Sesma et al. [2011a] to advance a new theory
for H/V. Based on the diffuse field assumption (DFA), Sánchez-Sesma et al. [2011a] linked
the H/V of the auto-correlated signal to the ratio of the imaginary parts of the Green’s
functions. The observed H/V can be compared to its theoretical counterpart, and provides
a basis for estimating the subsurface structure. It represents an opportunity to obtain the
local subsurface velocity model with only one short 3-component measurement of the
ambient field (e.g., only 20 minutes in Spica et al. [2015]). This capability is of particular interest in geotechnical engineering, seismic exploration, and engineering seismology. To date it has only been applied to limited data sets and geological settings [Salinas
et al., 2014; Spica et al., 2015; Lontsi et al., 2015; Piña-Flores, 2015; Rivet et al., 2015;
Lontsi et al., 2016; García-Jerez et al., 2016; Piña-Flores et al., 2016; Perton et al., 2017].
While these initial results are promising, there is still a need to explore the potential of the
method and to validate the results against independent information.
In this contribution, we pursue the idea of Lontsi et al. [2015] to use multiple H/V
measurements obtained from receivers at depth in a borehole to estimate complex geological structure at a site. The Groningen area, is a well-studied area and it has ∼70 borehole
sites, such that it represents an ideal natural laboratory to test the method.
We start with an overview of the geological horizons of interest that are sampled by
the measurements and some background on constraints on structure from previous studies (section 2). We then outline the theoretical background and seismic processing necessary to compute H/V (section 4), and the discrete wave number (DWN; [Bouchon, 1981])
method used to compute the H/V ratios for the forward problem in section 4.2. In section
5, we describe the estimation procedure for the velocity models at several borehole sites.
Finally, we discuss results and compare them with velocity models obtained independently
(section 6). In this paper we only present a few inversions because our focus is primarily
on implementation of the method and the validation of the results.
2 Geological Setting and Existing Velocity Model
The Groningen area is characterized by flat, low-lying topography with altitude close
to mean sea level. The horizons of interest (i.e., those sampled by our measurements)
are the Paleogene, Neogene, and younger deposits overlying the North Sea Supergroup
[De Mulder et al., 2003; Vos, 2015]. This ∼800 m thick layer of unconsolidated sediments
contains a large degree of vertical and lateral heterogeneity due to the influence of the last
three ice ages and associated sea level fluctuations. The deposits range from fluvial braid
plain sands to shallow marine (intertidal) and terrestrial deposits of soft clays to distinct
organic-rich peat formations. The uppermost sedimentary sequence is characterized by a
succession of fluvial, glacial, and marine deposits that are crosscut by deep subglacial features ("tunnel valleys"), which are filled with sands and clays and buried under younger
sediments. The North Sea Supergroup consists primarily of alternating marine grey sands,
sandstones, and clays.
Our study benefits from important, independent information in the form of an integrated VS model from the surface to the base of the North Sea Supergroup [Kruiver et al.,
2017]. This 3D VS model is a synthesis of three different VS models obtained through a
comprehensive set of geological, geotechnical, and geophysical observations; each of them
having a different depth range and spatial resolution. The shallowest part of the model extends from the surface to ∼50 m below sea level, and is based on a high resolution 3D
geological model GeoTOP [Stafleu et al., 2011; Stafleu and C, 2016], combined with VS
distributions with depth for the sediments [Kruiver et al., 2017]. At depths of ∼40 m to
∼120 m below sea level, the model is constrained through inversion of Rayleigh wave
observations from spatially extensive reflection seismic surveys. The deepest part of the
model is based on the Pre-Stack Depth Migration velocity model derived from sonic logs.
It ranges from ∼70 m below sea level to the base of the North Sea Supergroup at (∼800
m). These three distinct models were spliced into a single model that covers the full depth
–3–
Confidential manuscript to be submitted to JGI
A
G03
G08
G61
G18
G23
G29
G34
G40
G51
G56
B
NL
BE
Figure 1.
DE
A) KNMI-NAM permanent borehole sites (the G-array) in the Groningen gas field area. The
inverted colored triangles refer to the H/Vs shown in Fig 2. The blue contour depicts the outline of the gas
field. White areas depict urban centers. B) The location of the Groningen gas field (red box) in the northern
part of the Netherlands.
–4–
Confidential manuscript to be submitted to JGI
range for site response analyses. In what follows, we refer to this velocity model as the
Deltares-NAM model. Other shallow velocity models were provided by Noorlandt et al.
[In review].
We use 1D VS profiles of the Deltares-NAM model at the borehole sites to validate
our analysis. It is important to keep in mind, however, that the resolution and sensitivity
of the methods used to construct the Deltares-NAM model and the H/V velocity profiles
differ, as discussed in section 4.3.
3 Data
We used continuous data from the shallow borehole network installed in the Groningen area (Fig. 1; the G-array). Each borehole site includes one accelerograph at the surface that is collocated with 4 geophones at 50 m depth intervals (-50, -100, -150, -200
m). The continuous data are freely available on the Royal Dutch Meteorological Institute
(KNMI) website. We downsample the continuous data to 50 sps, remove instrument response, and convert the time series to velocity prior to processing.
4 H/V for a Diffuse Wavefield
By definition, the H/V spectral ratio corresponds to the square root of the spectral
energy ratio of the horizontal amplitudes (with indices 1, and 2) over the vertical direction
(index 3) [Arai and Tokimatsu, 2004] :
s
E1 (x,ω) + E2 (x,ω)
H
(x,ω) =
.
V
E3 (x,ω)
Perton et al. [2009] showed that at an observer location (surface or depth), these spectral
energies (i.e., the directional energy densities, as defined in Perton et al. [2009]) are proportional to the average auto-correlations of the diffuse wavefield components, which in
turn are proportional to the imaginary parts of the Green’s function (Im[G]) [SánchezSesma et al., 2011b]:
D
E
Ei (x,ω) = ρω2 ui (x,ω)ui∗ (x,ω) ∝ −ωIm Gii (x, x,ω) .
where ω is the angular frequency, ui (x,ω) is the displacement field in the i direction at a
point x, Im[ ] stands for the imaginary part, Gii (x, x,ω) is the displacement Green’s function in the direction i at a point x due to the application of a unit point force in the same
direction applied at the same point. The symbol ∗ stands for the complex conjugate operator and the product ui (x,ω)ui∗ (x,ω) corresponds in the frequency domain to the autocorrelation of the displacement field in the time domain. The brackets h i represent averaging
over time. In what follows, the dependence on ω and x is implicit. The ∝ means that the
expressions are proportional by a factor that is independent of ω and x.
Eq. 2 represents the same case as classic ambient seismic field cross-correlations
[e.g., Weaver and Lobkis, 2004] but for the special case when the source and receiver are
the same. Within this theoretical framework, Sánchez-Sesma et al. [2011a] proposed a theoretical description of H/V ratios and suggested that the H/V spectral ratio recorded at a
receiver could also be computed in terms of the imaginary part of the GF:
v
s
t
hu1 (x,ω) 2 i + hu2 (x,ω) 2 i
H
Im(G11 + G22 )
(x,ω) =
=
.
2
V
Im(G33 )
hu3 (x,ω) i
Eq. 3 links the average energy densities (i.e. the ambient field measurements; see section
4.1) with the Green’s function (i.e., the theoretical counterpart; see section 4.2) and treats
the H/V spectral ratio as an intrinsic property of the medium. It naturally allows for the
inversion (see section 5) of H/V that includes contributions of the full wavefield - i.e., including both surface and body waves.
–5–
(1)
(2)
(3)
Confidential manuscript to be submitted to JGI
4.1 Observed H/V
In the context of the DFA, eq. 2 is only valid when the seismic wave field is equipartitioned, i.e. all the incident waves are (e.g., P, S, or Rayleigh waves) have the same
energies [Perton et al., 2016]. This assumption is difficult to verify and unlikely to be true
in the majority of ambient field data. Therefore, some signal processing must be applied
to enhance the equipartitioning of the seismic wavefield, just as for traditional ambient
seismic field cross-correlation [e.g., Bensen et al., 2007]. As in Spica et al. [2015] and
Perton et al. [2017], we apply spectral whitening, which corresponds to source deconvolution. Because several sources can act in different frequency bands, the operation consists of normalizing the signals by the source energies computed
in each time window
q
P3
2
and across several frequency bands as ṽi (x,ω) = vi (x,ω)/
i=1 |vi (x, ∆ω)| ; where ∆ω
is a frequency band of 2.5 H z width centered on ω. We work with the particle velocity
v j (x,ω) = iωu j (x,ω) to preserve the link to energy. To remove only the spectral envelope,
the bandwidth has to be much larger than the oscillations in the spectra. The H/V spectral
ratio is therefore computed in terms of the wavefield auto-correlations as:
H
(x,ω) =
V
v
t
hv1 (x,ω) 2 i + hv2 (x,ω) 2 i
.
hv3 (x,ω) 2 i
Eq. 4 requires that the averaging is performed separately for each component. In that
sense, it is different from the calculation of the usual H/V spectral ratio [Nakamura, 1989],
which corresponds to the average of the ratios: H/V = hHw /Vw i. The average is computed on the spectra obtained over one day of continuous data that is windowed into sections of 100 s duration with an overlap of 20%. Each time window is demeaned, detrended, and bandpass filtered from 0.1 to 10 Hz.
The attenuation is generally a severe limitation for a field to become diffuse and
then equipartitioned. Several studies have shown this limitation for cross-correlation with
large inter-station distance [e.g., Lawrence and Prieto, 2011]. However, here, the interstation distance is null and the attenuation has a similar effect on waves in all directions.
Provided there is sufficient nearby sources of noise, attenuation should not be a limitation.
We have also shown theoretically that in presence of attenuation the auto-correlation is
proportional to the Im[G] for the same media but without attenuation [Perton and SánchezSesma, 2014]. A possible limitation of the method is, however, the presence of strong
anisotropy in the shallow sediment. This is not taken into account in this contribution.
Both the evaluation of anisotropy in the area and its assessment through H/V spectral ratio
are possible future research directions.
–6–
(4)
Confidential manuscript to be submitted to JGI
H/V
4
G03
G61
G18
G23
2
Depth [m]
04
100
H/V
G08
G29
100
G34
100
G40
100
G51
100
G56
2
0
100
freq. [Hz]
100
freq. [Hz]
100
freq. [Hz]
100
freq. [Hz]
100
freq. [Hz]
Figure 2.
H/V at different depth levels computed along the NW-SE line of stations shown as inverted col√
ored triangles in Fig. 1. The dashed gray lines in each panel is the 2 level (see section 4.3). The H/V in the
green frames (G03, G23, G29 and G56) are shown in Fig. 5 along with their inversion results.
Fig. 2 shows different H/V curves computed for the NW-SE line of boreholes highlighted in Fig. 1. By far, the strongest variations in the H/V curves are observed for the
surface measurements (black lines in Fig. 2). We attribute the higher amplitude of the
surface H/V ratios to a stronger impedance contrast the air-solid interface and a stronger
contribution from surface waves. The H/V curves at depth have much lower amplitude
variations, but slight changes in their shape reflect structural changes at depth. These are
expressions of shallow structural variability from north-west to the south-east in the study
area.
In order to test that the observed H/V is not biased due to lack of diffusivity of the
ambient field, we analyse the variability between H EV−W and H NV− S for stations at the surface (Fig. S.1). The fact that these measurements are similar at almost all the sites, suggests a near isotropic illumination by the noise sources [Perton et al., 2017]. Equipartition
of the wavefield is also discussed in section 4.3.
4.2 Theoretical H/V
In the second term of eq. 3, the Im[G] components must be related to a geometry and material properties that explain the data. Perton et al. [2017] explored the case of
strong lateral heterogeneities along a crater cliff; however, the case of Groningen does not
have strong topographic variation of the surface, which allows us to use a 3D unbounded
layered elastic media and smooth variation of the subsoil in the vicinity of the stations as
a reasonable starting assumption.
Several methods have been proposed to model the Gii components of H/V in a layered medium under the DFA [Sánchez-Sesma et al., 2011b; García-Jerez et al., 2016; Perton and Sánchez-Sesma, 2016; Lontsi et al., 2015]. Here, we use the Discrete Wave Number method (DWN; [Bouchon, 2003, 1981]. The DWN method is efficient and suitable for
solving the forward problem for H/V for stations located at the surface [Sánchez-Sesma
et al., 2011a; Perton et al., 2017]. In what follows, we discuss its efficiency when several
stations are located in a single geometrical and material configuration at different vertical
and/or horizontal positions.
As in [Sánchez-Sesma et al., 2011a], the wavefield is decomposed as a sum of plane
waves according to the horizontal component of the wavenumber vector in the radial direction k r . The Green’s function is the sum of two contributions, one due to the P and
–7–
0
50
100
150
200
Confidential manuscript to be submitted to JGI
SV waves and another due to SH waves. A compact form of the displacement P-SV Green’s
function component Giin (x) inside the layer n, where receiver and point source are superimposed at x = {x, z} and where the vector source is oriented in the same direction i is:
Giin (x) = GiiF Sn (x) +
Z
+∞
−∞
S %Pn ui%Pn + S %SV n ui%SV n + S &Pn ui&Pn + S &SV n ui&SV n dk r
(5)
For an extended description of the integrand, the interested reader can consult the textbook
of Aki and Richards [2002]. G F Sn is the incident field as in a full space (F S) having the
same properties as the layer n. The wave amplitudes S are unknowns that are solved for
at discrete values of k r by enforcing the continuity of displacement and traction at each
interface. Superscripts P and SV are for P- and SV-waves. We write the system of equations in matrix form [A][S] = [B]; where [S] = {..., S %Pn , S &Pn , S %SV n , S &SV n , ...} is the
unknown vector and [B] is the vector associated with the stress and displacement components of the source (i.e. to GiiF Sn ). [A] is comprised of the amplitude coefficients of the
boundary conditions and has dimensions equal to [4(N − 1)) ∗ (4(N − 1)] with N being the
number of layers.
As several sources are considered at different depths and for different directions, the
system of equations should be solved several times (i.e., for each source); however, since
we consider a single configuration (i.e., the matrix [A] is the same), it is more efficient to
calculate [A] −1 just once. Therefore, the wave amplitudes related to each source are calculated as the result of the multiplication [S] = [A] −1 [B], rather through inversion by Gaussian elimination for each vector [B]. In this way all the Green’s function’ components Giin
can be efficiently evaluated at multiple receiver positions and for the several directions.
The same procedure applies for the Green’s function associated with the SH contribution.
This allows consideration of all H/V for the same borehole site simultaneously during the
joint inversion.
–8–
Confidential manuscript to be submitted to JGI
4.3 Forward modeling using a complex velocity model
C1
C2
Gradient
Figure 3.
Theoretical H/V computed from the Deltares-NAM model (blue lines) and observed H/V (black
lines) at borehole G03. The letters C1 and C2 refer to possible strong velocity contrasts in the upper 50 m
depth range, and the dashed line refers to the gradient-like part of the velocity model.
We show in Fig. 3 the Deltares-NAM VS model at borehole site G03. In the 50 first
meters of the model, there are two strong velocity contrasts (C1 and C2 in Fig. 3). The
deeper velocity dependence is mainly expressed by a velocity gradient (dashed black line
in Fig. 3). Observed (black lines) and computed H/V (blue lines) at each depth level are
shown on the right part of the figure.
The discrepancies between the observed and synthetic H/V suggest that the inversion
of the observed H/V will result in a different velocity model, but with certain common
characteristics. While the model of Kruiver et al. [2017] is site-specific and has a resolution that changes with depth (see section 2), the H/V are sensitive to elastic properties
over a larger horizontal area that depends on the wavelength of the waves that contribute
to the observations [Perton et al., 2009; Piña-Flores et al., 2016]. The joint inversion result
will therefore give an averaged (vertical and horizontal) velocity model at a borehole site.
The H/V at the surface show one clear frequency peak, while H/V at depth shows
small oscillations superimposed on a flat trend. As recently detailed by Piña-Flores et al.
[2016], H/V is sensitive to both surface and body waves. The H/V spectra from receivers
at the free surface mainly originate from the strong contribution of surface waves and
show important localized frequency peaks of high amplitude. These frequency peaks are
dependent on the existence of strong impedance contrasts at depth. In the case of a simple
velocity model of one layer over a half-space (with shear velocity β1 and thickness h1 ),
the main peak frequency f bet a1 is well approximated by [Yamanaka et al., 1994; Piña-
–9–
Confidential manuscript to be submitted to JGI
Flores et al., 2016]:
β1
.
4h1
When two layers are considered on top of a half space, as suggested by the DeltaresNAM model (i.e., C1 and C2 in Fig. 3), H/V can present two peaks depending on the
thickness of each layer and the impedance contrasts at their boundaries. If the thickness
of the deepest layer is greater than the shallower layer and a sufficient impedance contrast
between these three materials exists, the two main peaks are well separated [BonnefoyClaudet et al., 2006; Field and Jacob, 1995; Bard and SESAME-team, 2004; Piña-Flores
et al., 2016]. Therefore, if β2 β1 and h1 h2 , the shear velocity of the superficial
layer can also be evaluated using eq. 6 (here f β1 is the high frequency peak) and the approximate frequency associated with the peak at low frequencies can be evaluated from:
f β1 =
f β12 =
1
4( hβ11
+
h2
β2 )
;
where β2 and h2 are the thickness and shear velocity of the second layer. While more
complex formulae that reflect several major effects of the model on the resonance frequency exist [Tuan et al., 2016], modeling of eq. 7 appears to be a good approximation
for layered structure [e.g., Piña-Flores et al., 2016]; however, since H/V at the free surface
in Fig. 3 does not show such character (i.e., a double frequency peak), we would not expect the inverted velocity model to develop two strong velocity contrasts in the first 50 m
at this site.
Because surface wave energy quickly decreases with depth, H/V in boreholes becomes more sensitive to body waves. While surface waves propagate in 2D space and
are generally not strongly reflected by lateral heterogeneity, body waves propagate in 3D
space and are reflected by the free surface and also by possible strong impedance contrasts at depth. As shown in [Perton et al., 2009] for a half space, the waves that travel
vertically up and down interfere and result in spectral oscillation periods in the energy
density components (E1 , E2 , E3 ). The amplitude of these oscillations decays with frequency and√ depth. In a half space, the dependence on frequency and depth are similar.
√
The H = E1 + E2 is mainly sensitive to the shear wave velocity while the V = E3 is
mainly sensitive to the compressional wave velocity. Also, as a requirement of equipartition, the three energy densities tend to the same value with depth.√Therefore, from eq. 3,
H/V in a half space should present these oscillations and tend to 2 at high frequency and
at depth (as observed in Fig. 2).
The gradient-like velocity structure suggested by the Deltares-NAM model does not
modify these conclusions since the observed H/V at depth expresses these features. The
latter suggests that we might be able to evaluate the velocity at each sensor depth; however, the absence of strong impedance contrast at depth reduces our ability to recover a
wave reflected away from the buried station and H/V is only weakly sensitive to structure
deeper than the deepest sensor (i.e., ≥200 m).
Other examples of observed H/V based on the complex velocity model provided by
Deltares-NAM are shown with the inversion results below. This allows us to assess how
close the site-specific model is from the H/V measurements and to assess the prospects for
improvement.
5 Inversion
As discussed in Piña-Flores et al. [2016], consideration of H/V at surface alone is
insufficient to characterize shallow properties uniquely since velocities and thicknesses
trade-off and lead to a similar H/V. Additionally, the forward problem is highly non-linear
and depends on several uncorrelated parameters [Piña-Flores et al., 2016; García-Jerez
et al., 2016]. Constraining the inversion by adding observations from sensors at depth significantly reduces the possible range of parameters [Lontsi et al., 2015]. Starting with an
accurate first guess - as for example a guess obtained from an independent measurement
–10–
(6)
(7)
Confidential manuscript to be submitted to JGI
such as a dispersion curve [Scherbaum et al., 2003; Piña-Flores, 2015; Lontsi et al., 2016]
- allows more rapid convergence. Although the Deltares-NAM model exists, our goal is to
provide independent measurements in our analysis to test the validity of the method.
5.1 Starting models from borehole interferometry
H1
Z=0
Z=50
Z=100
Z=150
Z=200
Figure 4.
H2
rs
rb
S
S
rb
S
rb
S
S
A) Geometry of the boreholes where the gray stars represent the virtual sources (S) and triangles
represent boreholes (r b ) or surface receivers (r s ). Ambient field cross-correlations are performed for every
depth level independently. B) Obtained correlation functions for the north-south and east-west components.
The picked travel time velocities (orange dots) are marked for every depth level and the mean velocity for the
entire well is marked in red. Note that similar travel times are picked in both horizontal components, suggesting weak horizontal anisotropy at this site. C) Resulting average velocity model (init) used as starting model
for inversion, compared to the Deltares-NAM velocity model harmonically averaged over 50 meter intervals.
We use ambient field borehole interferometry between pairs of adjacent overlying
sensors [e.g., Miyazawa et al., 2008] to obtain the mean shear wave travel time. In Fig.
4A, we show the geometry of the borehole stations and the sensor pairs used for interferometry. The cross-correlation of the ambient seismic field v(z S ) and v(zr ) in the frequency domain is given by:
C AB (ω) = hv(z S ) ∗ v(zr )i,
where z S and zr are the depths of the virtual source and of the receiver, respectively,
which are any of the consecutive sensors. Then, the imaginary part of the Green’s function between pairs of sensors (Im[G](z S , z r ) ) is obtained from the average correlation using:
hv(z S ) ∗ v(zr )i
Im[G](z S , z r ) ∝
;
|S(ω)| 2
where |S(ω)| 2 is the power spectrum of the noise [e.g., Wapenaar and Fokkema, 2006].
–11–
(8)
(9)
Confidential manuscript to be submitted to JGI
The cross-correlations are computed for both east-west and north-south component
pairs between the wave motion observed at a depth z S and another directly overlying sensor zr . Only ambient noise is used and small earthquakes are removed [e.g., Miyazawa
et al., 2008] based on the KNMI earthquake catalogue. Correlations are computed over
one month of continuous data for 25 s windows. Fig. 4b shows the resulting bandpass filtered (1.5-8 Hz) stacked correlations in the time domain. An impulsive arrival is observed
on both causal and anti-causal parts of the Green’s function, corresponding to upward and
downward propagating S-waves.
We estimate the arrival time as in Nakata and Snieder [2012]; i.e., by seeking the
three adjacent samples with the largest amplitude values and applying a quadratic interpolation to find the time at which the wave has maximum amplitude. This time is taken
as the travel time for a shear wave that propagates between the borehole and overlying
sensor. We use the average of estimated travel time on the east-west and north-south components to estimate an average velocity layer. The obtained velocity model at site G43 is
shown in Fig. 4c and is compared to the Deltares-NAM velocity model harmonically averaged over 50 meter intervals.
These results might be biased by multiple arrivals caused by shear wave reflection
inside the 50 m section due to strong geological contrasts. Besides, if the noise is not
fully diffuse, Green’s functions are not well retrieved and some small bias can be observed
in the travel times [Tsai, 2009]. In such case, the interferometry picked travel time would
be smaller than the expected one, which results in higher velocities. Because the average
velocity models tend to underestimate the Deltares-NAM velocity model at most of the
sites (i.e., Fig. 4c and Fig. S.2), and because the analysis performed on the H/Vs in Fig.
S.1, this is not likely to be a strong effect. The discrepancy observed in the Fig. S.2 is
also partly explained by the fact that the interferometry picked travel-times are directly
translated to average velocity. Contrary to the travel-time, the velocity estimated by the
interferometry should not be exactly equal to the average of the velocities in the 50 m
section since it corresponds to a harmonic average.
5.2 Parametrization, misfit function and inversion
The soil properties and the layer thicknesses inside a borehole site are assessed by
inverting jointly H/V at several depths, z. The objective function is defined as the root
mean square difference between observed and predicted H/V computed with the DWN
method:
v
u
t
!2
1 X X H e x p H DW N
−
;
=
Nz Nω z ω V
V
i
where Nz = 5 is the number of depths at which z i considered in the joint inversion, and
Nω = 14 is the number of points taken in the spectra from 1. to 8. H z. Fluctuations in the
H/V spectra at higher frequency are associated with small, local heterogeneities, which are
of minor importance here [Piña-Flores et al., 2016].
The only free parameters considered during the inversion are the shear wave velocities ( β). To reduce the number of parameters, the compressional velocity (α) and the mass
density are assumed to be related to the shear velocity through polynomial relationships
[Brocher, 2005; Berteussen, 1977]. We impose that one velocity layer in each section of
50 m (i.e. between each sensor pair) is constrained by the interferometry picked travel
time and by the other velocities of the section. In other words, if i is the layer index of
height hi and shear velocity βi , and j0 ≤ j ≤ j1 are the indices of all the layers belonging
to the 50 m section number k then βi should be comprised between [0.85 − 1.15]hi /t i
P
with t i = t kcor r − j={ j0 : j1 }, j,i h ∗j / β j . For j0 and j1 , only part of the thicknesses that belongs to the section is considered. Because this constraint strongly relies on assignment
of the average velocity from interferometry, we allow the constrained velocity to vary by
15% during the inversion. Also, because of this constraint, the inverted velocity model
will always be closer to the starting velocity model than to the Deltares-NAM model. The
thicknesses are assumed constant during the inversion, and the model is refined iteratively
–12–
(10)
Confidential manuscript to be submitted to JGI
[e.g., Spica et al., 2016, 2017]. The layers that are refined between two iterations are chosen based on the sensitivity of the misfit to a small velocity change .
We use a Pattern-Search method for the iterating inversion because of the strong
non-linearity of the problem and because of its efficiency [e.g., Audet and Dennis Jr, 2002].
All the inversions presented in this paper have a misfit value () lower than 0.01.
6 Results and discussion
We show examples of joint inversion results at borehole sites G03, G18, G34, and
G56 in Fig. 5. The left panels depict observed H/V (black lines) along with the best fits
after inversion (green lines) and H/V computed from the complex Deltares-NAM velocity
model (blue lines). The right panel depicts the best velocity model (green line) along with
initial model obtained from borehole interferometry (dashed black line) and the DeltaresNAM model (blue lines). Results at the other borehole sites are shown in the auxiliary
material (Fig. S.3).
The inverted velocity models are in good agreement with the results of [Kruiver
et al., 2017]. For the first 50 m section, the inverted velocity models at G03, G18 and
G34 show the same trends as the Deltares-NAM velocity model, i.e. some high contrast
with a higher local velocity at approximately the same depth; however, the velocity ratio
between these two models is sometimes close to 2. This difference could be explained by
the use of the relationship between β, α and ρ, which could be incorrect at shallow depth
where the soil is likely fluid-saturated and where shear velocity might be expected to drop
dramatically. Allowing the compressional velocity to become a free parameter, however,
would result in an unmanageable increase in the number of unknown parameters for the
inversion.
Below 100 m depth, the Deltares-NAM and inverted H/V models agree well, and
they are concordant with the gradient-like trend. This good agreement suggests that the
different volumes sampled by the two studies have similar properties. The use of H/V at
depth in the inversion problem appears to reduce the non-uniqueness significantly and, potentially, to successfully reconstruc a accurate velocity model. As mentioned in section
4.3, some of the discrepancies observed between the different sets of velocity models may
be attributable to the different sensitivities of the two approaches. An important difference is the greater number, and strength of velocity reductions with increasing depth. This
could be due to the strong constraints imposed by the interferometry picked travel time
during the inversion. On the other hand, we note that the theoretical H/V computed from
the Deltares-NAM model fits the observed H/V fairly well, especially at shallow depths.
This may suggest that we are converging to an accurate velocity model and hence an accurate site characterization.
–13–
Confidential manuscript to be submitted to JGI
Figure 5.
Inversion results at sites G03, G18, G34 and G56. The left panels of each subplot depict the
observed H/Vs (black lines) along with the inverted H/Vs (green lines), the theoretical H/Vs corresponding
to the complex Deltares-NAM velocity model (blue lines) and the H/V spectral ratio uncertainty range (pink
dashed lines). The upper and lower bounds of the auto-correlations are represented by black dashed lines.
They are obtained as the maximum and minimum at each frequency from all the auto-correlations computed
with 50 windows numbers, which allows convergence. The related velocity models are shown in the right
panels of each subplot.
–14–
Confidential manuscript to be submitted to JGI
7 Conclusions
We present new theoretical and empirical results on the computation and the inversion of H/V spectral ratios based on the diffuse field assumption for combinations of
receivers at surface and depth. First, we obtained the mean shear travel times by applying ambient field interferometry between adjacent sensors inside the boreholes. By crosscorrelating the ambient seismic field we were able to retrieve upward and downward propagating S-waves between adjacent borehole seismic stations. We used this mean velocity model as a starting velocity model for the joint inversion of the depth-dependent H/V.
We found it especially useful to constrain the inversion for velocities deeper than 100 m.
The use of 5 independent H/V measurements at different depths, together with the average
shear wave travel times, helps reduce the range of acceptable parameters and helps speed
convergence of the solution. In this inversion, the theoretical H/V counterparts (i.e., the
forward problem) were computed using the DWN method, which allowed efficient modeling of the wave field for different sources at the surface and at depth.
We successfully obtained complex VS velocity profiles of the shallow sub-surface
at different borehole sites in the Groningen area. Velocity models are globally in good
agreement with previous site characterization for the region [Kruiver et al., 2017]. Our
approach has the potential to reduce uncertainty in modeling the response of the shallow
crust, which is an important component of probabilistic seismic hazard analysis at Groningen. Our results also validate the 3D character of the ambient field sources and motivate
measurements of H/V at depth. They also demonstrate the power and reliability of the
method, which could be applied elsewhere, where other constraints on shallow structure
are lacking. The methodology presented here allows reliable recovery of layered velocity
structure at boreholes sites.
Acknowledgments
We would like to thank Pauline Kruiver from Deltares for providing us with their
velocity models at borehole sites [Kruiver et al., 2017]. We thank the KNMI to make
data accessible. We appreciate the contribution of two anonymous reviewers. All the figures have been plotted with matplotlib [Hunter, 2007] and some of the data processing
steps have been performed using obspy [Beyreuther et al., 2010] and pyrocko (available at:
http://pyrocko.org/).
–15–
Confidential manuscript to be submitted to JGI
H/V
3
2
G08
G61
G18
G23
1
0
3
H/V
G03
2
G29
G34
G40
G51
G56
1
0
1 3 5 7
Figure S .1.
1 3 5 7
H E −W
V
1 3 5 7
frequency [Hz]
1 3 5 7
1 3 5 7
and H NV− S for stations at the surface. Same stations as in Fig. 1 and 2.
–16–
HN − S
V
HE − W
V
Confidential manuscript to be submitted to JGI
Figure S .2.
Comparison between the Deltares-NAM velocity model and the initial model form interferom-
etry. Most of the initial models have lower velocities in the first 100 m.
–17–
Confidential manuscript to be submitted to JGI
Figure S .3. Inversion results at all borehole sites highlighted in Fig. 1.
–18–
Confidential manuscript to be submitted to JGI
References
Aki, K. (1957), Space and time spectra of stationary stochastic waves, with spectral reference to microtremors, Bull. Earthq. Res. Inst., 35, 415–456.
Aki, K., and P. Richards (2002), Quantitative seismology, University Science Book.
Arai, H., and K. Tokimatsu (2004), S-wave velocity profiling by inversion of microtremor
H/V spectrum, Bulletin of the Seismological Society of America, 94(1), 53–63, doi:10.
1785/0120030028.
Audet, C., and J. E. Dennis Jr (2002), Analysis of generalized pattern searches, SIAM
Journal on Optimization, 13(3), 889–903.
Bard, P.-Y., and SESAME-team (2004), Guidelines for the implementation of the H/V
spectral ratio technique on ambient vibrations measurements, processing and interpretations, SESAME european research project EVG1-CT-2000–00026, deliverable d23.12,
available at http://sesame-fp5.obs.ujf-grenoble.fr.
Bensen, G. D., M. H. Ritzwoller, M. P. Barmin, A. L. Levshin, F. Lin, M. P. Moschetti,
N. M. Shapiro, and Y. Yang (2007), Processing seismic ambient noise data to obtain
reliable broad-band surface wave dispersion measurements, Geophysical Journal International, 169(3), 1239–1260, doi:10.1111/j.1365-246X.2007.03374.x.
Berteussen, K. A. (1977), Moho depth determinations based on spectral-ratio analysis of
NORSAR long-period P waves, Physics of the Earth and Planetary Interiors, 15(1), 13–
27, doi:10.1016/0031-9201(77)90006-1.
Beyreuther, M., R. Barsch, L. Krischer, T. Megies, Y. Behr, and J. Wassermann (2010),
ObsPy: A Python Toolbox for Seismology, Seismological Research Letters, 81(3), 530–
533, doi:10.1785/gssrl.81.3.530.
Bommer, J. J., P. J. Stafford, B. Edwards, B. Dost, E. van Dedem, A. Rodriguez-Marek,
P. Kruiver, J. van Elk, D. Doornhof, and M. Ntinalexis (2017), Framework for a
ground-motion model for induced seismic hazard and risk analysis in the Groningen
gas field, the Netherlands, Earthquake Spectra.
Bonnefoy-Claudet, S., C. Cornou, P.-Y. Bard, F. Cotton, P. Moczo, J. Kristek, and D. Fah
(2006), H/V ratio: a tool for site effects evaluation. results from 1-D noise simulations,
Geophysical Journal International, 167(2), 827–837, doi:10.1111/j.1365-246X.2006.
03154.x.
Bonnefoy-Claudet, S., A. Köhler, C. Cornou, M. Wathelet, and P.-Y. Bard (2008), Effects
of Love waves on microtremor H/V ratio, Bulletin of the Seismological Society of America, 98(1), 288–300, doi:10.1785/0120070063.
Bouchon, M. (1981), A simple method to calculate Green’s functions for elastic layered
media, Bulletin of the Seismological Society of America, 71(4), 959–971.
Bouchon, M. (2003), A review of the discrete wavenumber method, Pure and applied Geophysics, 160(3-4), 445–465.
Brocher, T. M. (2005), Empirical relations between elastic wavespeeds and density in the
Earth’s crust, Bulletin of the Seismological Society of America, 95(6), 2081–2092.
Capon, J. (1969), High-resolution frequency-wavenumber spectrum analysis, Proceedings
of the IEEE, 57(8), 1408–1418, doi:10.1109/PROC.1969.7278.
De Mulder, E. F., M. C. Geluk, I. Ritsema, W. E. Westerhoff, and T. E. Wong (2003), De
ondergrond van Nederland, Geologie van Nederland, deel, 7.
Field, E. H., and K. H. Jacob (1995), A comparison and test of various site-response estimation techniques, including three that are not reference-site dependent, Bulletin of the
Seismological Society of America, 85(4), 1127–1143.
García-Jerez, A., J. Piña-Flores, F. J. Sánchez-Sesma, F. Luzón, and M. Perton (2016), A
computer code for forward calculation and inversion of the H/V spectral ratio under the
diffuse field assumption, Computers & Geosciences, 97, 67 – 78, doi:http://dx.doi.org/
10.1016/j.cageo.2016.06.016.
Hunter, J. D. (2007), Matplotlib: A 2D graphics environment, Computing in Science and
Engineering, 9(3), 90–95, doi:10.1109/MCSE.2007.55.
–19–
Confidential manuscript to be submitted to JGI
Kruiver, P. P., E. van Dedem, R. Romijn, G. de Lange, M. Korff, J. Stafleu, J. L. Gunnink,
A. Rodriguez-Marek, J. J. Bommer, J. van Elk, et al. (2017), An integrated shear-wave
velocity model for the Groningen gas field, the Netherlands, Bulletin of Earthquake Engineering, pp. 1–26.
Lacoss, R., E. Kelly, and M. Toksöz (1969), Estimation of seismic noise structure using
arrays, GEOPHYSICS, 34(1), 21–38, doi:10.1190/1.1439995.
Lawrence, J. F., and G. A. Prieto (2011), Attenuation tomography of the western United
States from ambient seismic noise, J. Geophys. Res, 116, B06,302.
Lontsi, A. M., F. J. Sánchez-Sesma, J. C. Molina-Villegas, M. Ohrnberger, and F. Krüger
(2015), Full microtremor H/V (z, f) inversion for shallow subsurface characterization,
Geophysical Journal International, 202(1), 298–312.
Lontsi, A. M., M. Ohrnberger, F. Krüger, and F. J. Sánchez-Sesma (2016), Combining
surface-wave phase-velocity dispersion curves and full microtremor horizontal-tovertical spectral ratio for subsurface sedimentary site characterization, Interpretation,
4(4), SQ41–SQ49.
Malischewsky, P. G., and F. Scherbaum (2004), Love’s formula and H/V-ratio (ellipticity)
of Rayleigh waves, Wave Motion, 40(1), 57–67, doi:10.1016/j.wavemoti.2003.12.015.
Miyazawa, M., R. Snieder, and A. Venkataraman (2008), Application of seismic interferometry to extract P-and S-wave propagation and observation of shear-wave splitting
from noise data at Cold lake, Alberta, Canada, Geophysics, 73(4), D35–D40.
Nakamura, Y. (1989), A method for dynamic characteristics estimation of subsurface using microtremor on the ground surface, Railway Technical Research Institute, Quarterly
Reports, 30(1).
Nakamura, Y. (2000), Clear identification of fundamental idea of Nakamura’s technique
and its applications., in Proceedings of the 12th world conference on earthquake engineering., New Zealand: Auckland.
Nakata, N., and R. Snieder (2012), Estimating near-surface shear wave velocities in Japan
by applying seismic interferometry to kik-net data, Journal of Geophysical Research:
Solid Earth, 117(B1).
Nakata, N., J. P. Chang, J. F. Lawrence, and P. Boué (2015), Body wave extraction and
tomography at Long Beach, California, with ambient-noise interferometry, Journal of
Geophysical Research: Solid Earth, 120(2), 1159–1173, doi:10.1002/2015JB011870.
Noorlandt, R. P., P. P. Kruiver, M. P. E. Kleine, G. Karaoulis, G. de Lange, A. Di Matteo,
J. von Ketelhodt, E. Ruigrok, B. Edwards, A. Rodriguez-Marek, J. Bommer, J. van Elk,
and D. Doornhof (In review), Characterisation of ground-motion recording stations in
the Groningen gas field, Journal of Seismology.
Olsen, K. (2000), Site amplification in the Los Angeles basin from three-dimensional
modeling of ground motion, Bulletin of the Seismological Society of America, 90(6B),
S77–S94.
Perton, M., and F. J. Sánchez-Sesma (2014), Normalization during the processing of noise
correlations, in IUGG, Merida.
Perton, M., and F. J. Sánchez-Sesma (2016), Green’s function calculation from equipartition theorem, The Journal of the Acoustical Society of America, 140(2), 1309–1318.
Perton, M., F. J. Sánchez-Sesma, A. Rodríguez-Castellanos, M. Campillo, and R. L.
Weaver (2009), Two perspectives on equipartition in diffuse elastic fields in three dimensions, The Journal of the Acoustical Society of America, 126(3), 1125–1130, doi:
10.1121/1.3177262.
Perton, M., M. A. Contreras-Zazueta, and F. J. Sánchez-Sesma (2016), Indirect boundary
element method to simulate elastic wave propagation in piecewise irregular and flat regions, Geophysical Journal International, 205(3), 1832–1842.
Perton, M., Z. Spica, and C. Caudron (2017), Inversion of the horizontal to vertical spectral ratio in presence of strong lateral heterogeneity, Geophysical Journal International,
under review.
–20–
Confidential manuscript to be submitted to JGI
Piña-Flores, J. (2015), Cálculo e inversión del cociente H/V a partir de ruido ambiental,
MSc thesis, Universidad Nacional Autónoma de México, Mexico city.
Piña-Flores, J., M. Perton, A. García-Jerez, E. Carmona, F. Luzón, J. C. Molina-Villegas,
and F. J. Sánchez-Sesma (2016), The inversion of spectral ratio H/V in a layered system
using the diffuse field assumption (dfa), Geophysical Journal International, doi:10.1093/
gji/ggw416.
Rivet, D., M. Campillo, F. Sanchez-Sesma, N. M. Shapiro, and S. K. Singh (2015), Identification of surface wave higher modes using a methodology based on seismic noise and
coda waves, Geophysical Journal International, 203(2), 856–868.
Roux, P., L. Moreau, A. Lecointre, G. Hillers, M. Campillo, Y. Ben-Zion, D. Zigone, and
F. Vernon (2016), A methodological approach towards high-resolution surface wave
imaging of the san Jacinto fault zone using ambient-noise recordings at a spatially
dense array, Geophysical Journal International, 206(2), 980–992.
Salinas, V., F. Luzón, A. García-Jerez, F. Sánchez-Sesma, H. Kawase, S. Matsushima,
M. Suarez, A. Cuellar, and M. Campillo (2014), Using diffuse field theory to interpret
the H/V spectral ratio from earthquake records in Cibeles seismic station, Mexico city,
Bulletin of the Seismological Society of America, 104(2), 995–1001.
Sánchez-Sesma, F. J., M. Rodríguez, U. Iturrarán-Viveros, F. Luzón, M. Campillo,
L. Margerin, A. García-Jerez, M. Suarez, M. A. Santoyo, and A. Rodríguez-Castellanos
(2011a), A theory for microtremor H/V spectral ratio: application for a layered medium,
Geophysical Journal International, 186(1), 221–225, doi:10.1111/j.1365-246X.2011.
05064.x.
Sánchez-Sesma, F. J., R. L. Weaver, H. Kawase, S. Matsushima, F. Luzón, and
M. Campillo (2011b), Energy partitions among elastic waves for dynamic surface loads
in a semi-infinite solid, Bulletin of the Seismological Society of America, 101(4), 1704–
1709, doi:10.1785/0120100196.
Scherbaum, F., K.-G. Hinzen, and M. Ohrnberger (2003), Determination of shallow shear
wave velocity profiles in the Cologne, Germany area using ambient vibrations, Geophysical Journal International, 152(3), 597–612.
Spica, Z., C. Caudron, M. Perton, T. Lecocq, T. Camelbeeck, D. Legrand, J. Piña-Flores,
A. Iglesias, and D. K. Syahbana (2015), Velocity models and site effects at Kawah Ijen
volcano and Ijen caldera (Indonesia) determined from ambient noise cross-correlations
and directional energy density spectral ratios, Journal of Volcanology and Geothermal
Research, 302, 173–189.
Spica, Z., M. Perton, M. Calò, D. Legrand, F. Córdoba-Montiel, and A. Iglesias (2016),
3-D shear wave velocity model of Mexico and south US: bridging seismic networks
with ambient noise cross-correlations (C 1 ) and correlation of coda of correlations (C 3 ),
Geophysical Journal International, 206(3), 1795–1813, doi:10.1093/gji/ggw240.
Spica, Z., M. Perton, and D. Legrand (2017), Anatomy of the Colima volcano magmatic
system, Mexico, Earth and Planetary Science Letters, 459, 1–13.
Stafleu, J., and D. W. C (2016), Product specification subsurface model geotop, version 1.3, Tech. Rep. R10133, 53 pp, TNO report, doi:https://www.dinoloket.nl/en/
subsurface-models.
Stafleu, J., D. Maljers, J. Gunnink, A. Menkovic, and F. Busschers (2011), 3D modelling
of the shallow subsurface of Zeeland, the Netherlands, Netherlands Journal of Geosciences, 90(04), 293–310.
Tsai, V. C. (2009), On establishing the accuracy of noise tomography travel-time measurements in a realistic medium, Geophysical Journal International, 178(3), 1555–1564,
doi:10.1111/j.1365-246X.2009.04239.x.
Tuan, T. T., P. C. Vinh, M. Ohrnberger, P. Malischewsky, and A. Aoudia (2016), An improved formula of fundamental resonance frequency of a layered half-space model used
in H/V ratio technique, Pure and Applied Geophysics, 173(8), 2803–2812.
Vos, P. C. (2015), Origin of the Dutch coastal landscape. Long-term landscape evolution of
the Netherlands during the Holocene, described and visualized in national, regional and
–21–
Confidential manuscript to be submitted to JGI
local palaeogeographical map series, Utrecht University.
Wapenaar, K., and J. Fokkema (2006), Green’s function representations for seismic interferometry, Geophysics, 71(4), SI33–SI46.
Weaver, R. L. (2005), Information from seismic noise, Science, 307(5715), 1568–1569,
doi:10.1126/science.1109834.
Weaver, R. L., and O. I. Lobkis (2004), Diffuse fields in open systems and the emergence
of the Green’s function (l), The Journal of the Acoustical Society of America, 116(5),
2731–2734, doi:http://dx.doi.org/10.1121/1.1810232.
Yamanaka, H., M. Takemura, H. Ishida, and M. Niwa (1994), Characteristics of longperiod microtremors and their applicability in exploration of deep sedimentary layers,
Bulletin of the Seismological Society of America, 84(6), 1831–1841.
–22–