A generalized theory for full microtremor horizontal-to-vertical [H/V (z, f)] spectral ratio interpretation in offshore and onshore environments
Agostiny Marrios Lontsi, Antonio García-Jerez, Juan Camilo Molina-Villegas, Francisco José Sánchez-Sesma, Christian Molkenthin, Matthias Ohrnberger, Frank Krüger, Rongjiang Wang, Donat Fäh
ppublished in
Geophys. J. Int.
A generalized theory for full microtremorhorizontal-to-vertical [ H/V ( z, f )] spectral ratiointerpretation in offshore and onshore environments Agostiny Marrios Lontsi (cid:63) , Antonio Garc´ıa-Jerez , , Juan Camilo Molina-Villegas , , Francisco Jos´e S´anchez-Sesma ,Christian Molkenthin , Matthias Ohrnberger , Frank Kr¨uger , Rongjiang Wang , Donat F¨ah Swiss Seismological Service, ETH Z¨urich, Switzerland, Departamento de Qu´ımica y F´ısica, Universidad de Almeria, Espa˜na, Facultad de Ingenier´ıas, Universidad de Medell´ın, Colombia, Instituto de Ingenier´ıa, Universidad Nacional Aut´onoma de M´exico, M´exico, Institute of Mathematics, University of Potsdam, Germany, Institute of earth and environmental Science, University of Potsdam, Germany, GFZ German Research Centre for Geosciences, Potsdam, Germany Instituto Andaluz de Geof´ısica, Universidad de Granada, Espa˜na, Departamento de ingenier´ıa civil, Facultad de Minas, Universidad Nacional de Colombia – Sede Medell´ın, Colombia.
Abstract
Advances in the field of seismic interferometry have provided a basic theoretical interpretation to thefull spectrum of the microtremor horizontal-to-vertical spectral ratio [ H/V ( f )] . The interpretationhas been applied to ambient seismic noise data recorded both at the surface and at depth. The newalgorithm, based on the diffuse wavefield assumption, has been used in inversion schemes to estimateseismic wave velocity profiles that are useful input information for engineering and exploration seis-mology both for earthquake hazard estimation and to characterize surficial sediments. However, until a r X i v : . [ phy s i c s . g e o - ph ] J u l A. M. Lontsi et al. now, the developed algorithms are only suitable for on land environments with no offshore considera-tion. Here, the microtremor
H/V ( z, f ) modeling is extended for applications to marine sedimentaryenvironments for a 1D layered medium. The layer propagator matrix formulation is used for the com-putation of the required Green’s functions. Therefore, in the presence of a water layer on top, thepropagator matrix for the uppermost layer is defined to account for the properties of the water column.As an application example we analyze eight simple canonical layered earth models. Frequencies rang-ing from . to Hz are considered as they cover a broad wavelength interval and aid in practice toinvestigate subsurface structures in the depth range from a few meters to a few hundreds of meters.Results show a marginal variation of percent at most for the fundamental frequency when a waterlayer is present. The water layer leads to variations in H/V peak amplitude of up to percent atopthe solid layers. Over the past decades, using the single-station microtremor horizontal-to-vertical (
H/V ) spectral ra-tio as a method for shallow subsurface characterization has attracted a number of site investigationstudies both on land (e.g. Bard, 1998; F¨ah et al., 2003; Scherbaum et al., 2003; Lontsi et al., 2015,2016; Garc´ıa-Jerez et al., 2016; Pi˜na-Flores et al., 2017; Spica et al., 2018; Garc´ıa-Jerez et al., 2019)and in marine environment (e.g. Huerta-Lopez et al., 2003; Muyzert, 2007; Overduin et al., 2015).The interest in the method is mainly due to its practicability, its cost efficiency, and the minimum in-vestment effort during microtremor (ambient noise or passive seismic) survey campaigns. The genericengineering parameter directly estimated from the spectrum of the microtremor
H/V spectral ratio isthe site fundamental frequency (e.g. Nakamura, 1989; Lachet & Bard, 1994). The fundamental fre-quency of a site generally corresponds to the frequency for which the microtremor
H/V spectral ratioreaches its maximum amplitude.Although the peak frequency is relatively well understood, this is not straightforward for sec-ondary peaks as they could represent higher modes or materialize the presence of more than onestrong contrast in the subsurface lithology. It is therefore important in the analysis to use a physi-cal formulation for the
H/V spectral ratio that not only accounts for the full spectrum (includingfirst and subsequent secondary peaks), but also includes all wave constituent parts. Based upon theadvances in seismic noise interferometry (e.g. Lobkis & Weaver, 2001; Shapiro & Campillo, 2004;Curtis et al., 2006; Wapenaar & Fokkema, 2006; Sens-Sch¨onfelder & Wegler, 2006; Gou´edard et al.,2008), S´anchez-Sesma et al. (2011) proposed a physical model for the interpretation of the full spec-trum of the microtremor
H/V spectral ratio. This has been extended to include receivers at depths(Lontsi et al., 2015). This additional information from receivers at depth is an added value during the icrotremor H/V(z, f) modeling in marine environment
H/V spectral ratio acquisition on land, no significant effort has been made forthe marine acquisition counterpart. An early study for a station on the seafloor was performed byHuerta-Lopez et al. (2003), assuming that the wavefield is due to the propagation of an incident planeSH body wave. With the evolving technology in borehole acquisition seismic instruments and datatransmission (e.g. Stephen et al., 1994), there is a growing need for efficient subsea exploration andgeohazard estimation as reported by Djikpesse et al. (2013).Here we further extend the diffuse field model (S´anchez-Sesma et al., 2011; Lontsi et al., 2015)to allow for the interpretation of the
H/V ( z, f ) both in marine sedimentary environment and on landeven though applicability to marine environments is emphasized.The Thomson-Haskell propagator matrix (Thomson, 1950; Haskell, 1953) is used to relate thedisplacement and stress for SH and P-SV waves at two points within an elastic 1D layered medium.The use of the propagator matrix formulation allows us to easily include a propagator for a layer on topthat accounts for the properties of the water layer and to subsequently compute the Green’s functionfor points at different depths. The classical Thomson-Haskell method is unstable when waves becomeevanescent. To remedy this issue, many attempts have been made (e.g Knopoff, 1964; Dunkin, 1965;Abo-Zena, 1979; Kennett & Kerry, 1979; Harvey, 1981; Wang, 1999). Here, we use the orthonor-malization approach by Wang (1999) which preserves the original Thomson-Haskell matrix algorithmand avoid the loss of precision by inserting an additional procedure that makes in-situ base vectorsorthonormal.A synthetic analysis is performed on eight simple canonical earth models. The models differ bythe presence of soft sediment structures with different overall thickness (two in total) and the presenceof a water column with varying depth at the top. The first sediment structure is a very simple onelayer over a half-space earth model and the second is a realistic structural model obtained from sitecharacterization at Baar, a municipality in the Canton of Zug, Switzerland. The H/V spectral ratiois estimated for frequencies ranging from . to Hz. The effects of the water column on the
H/V spectrum at selected depths are interpreted.
H/V
SPECTRAL RATIO: A PHYSICAL INTERPRETATION
Here, the main steps linking the microtremor
H/V ( z, f ) spectral ratio to the elastodynamic Green’sfunctions are presented. The basic expressions for SH and P-SV wave contributions to the Green’sfunctions and some considerations for numerical integration are summarized. A. M. Lontsi et al.
H/V ( z, f ) interpretation: Onshore case Starting from a three-component ambient vibration data, the microtremor
H/V spectral ratio at a givenpoint at the earth surface or at depth (onshore: Figure 1 without water layer) for a known frequency f is estimated using Equation 1. H/V ( z, f ) = (cid:115) E ( z, f ) + E ( z, f ) E ( z, f ) , (1)where E m ( z, f ) = ρω (cid:104) u m ( z, f ) u m ∗ ( z, f ) (cid:105) is physically regarded as the directional energy density, ρ is the mass density, ω is the angular frequency and u m ( m = 1 , , ) is the recorded displacementwavefield in the orthogonal direction m . The indexes m = 1 , correspond to the horizontal compo-nents while m = 3 corresponds to the index for the vertical component. The summation convention forrepeated indexes is not applied here. The symbol ∗ stands for complex conjugate. Using interferomet-ric principles under the diffuse field assumption, it can be shown that the average of the autocorrelationof the displacement field is proportional to the imaginary part of the Green’s function assuming thesource and the receivers are at the same point (S´anchez-Sesma et al., 2008; Snieder et al., 2009, see asummary in Appendix A). Equation 1 in terms of the Green’s function is expressed as: H/V ( z, f ) = (cid:115) Im [ G ( z, z, f )] + Im [ G ( z, z, f )] Im [ G ( z, z, f )] = (cid:115) Im [ G ( z, z, f )] Im [ G ( z, z, f )] (2)We are therefore left with the computation of the Green’s functions G = G and G . The elas-todynamic Green’s function in a 1D elastic layered medium (onshore: Figure 1 without water layer)is the set of responses for unit harmonic loads in the three directions. Using cylindrical coordinatesthe contribution of the radial-vertical (P-SV) and transverse (SH) motions are decoupled. Therefore,it suffices to solve each case separately using the integration on the horizontal wavenumber (Bouchon& Aki, 1977). Assuming the subsurface structure can be approximated by a stack of homogeneous layers over a half-space as depicted in Figure 1 where for example the j th layer is characterized in the onshore case bythe compressional wave velocity V P j , the shear wave velocity V S j , the density ρ j , the layer thickness h j , and the attenuation parameters Q P j and Q S j for the P- and and S-wave respectively, Im [ G ] ,Im [ G ] and Im [ G ] are given by: icrotremor H/V(z, f) modeling in marine environment . . .. . .. . . . . . Re ce iv e r p oin t a t d e p th . Wa te r1 V P1 , V S1 , ρ , h . . .. . .. . . . . . Wa te r . Re ce iv e r p o in t a t s u rfa ce V P2 , V S2 , ρ , h V Pj , V Sj , ρ j , h j V Pn , V Sn , ρ n , h n V P ∞ , V S ∞ , ρ ∞ , h ∞ V Pj , V Sj , ρ j , h j V Pn ,V Sn , ρ n , h n V P1 , V S1 , ρ , h V P2 , V S2 , ρ , h V P ∞ , V S ∞ , ρ ∞ , h ∞ h23jj+1nn+1 Figure 1.
Schematic representation of a 1D layered medium. The representation without the water layer on topcorresponds to the onshore case and the representation with water layer corresponds to the offshore case. For therepresentation on the left, the receiver location is at the earth surface when no water layer is present (onshore)and at the water (lake, sea, ocean) bottom when the water layer is present (offshore). For the representation onthe right, the receiver location is at depth. Except for the water layer in the offshore case where the shear wavevelocity is zero, any other layer j either onshore or offshore is characterized by the seismic parameters V P j , V S j , ρ j , h j , Q P j , and Q S j . Im [ G ] = 14 π (cid:90) ∞ Im [ g SH ] kdk + 14 π (cid:90) ∞ Im [ g PSV ] kdk (3)Im [ G ] = 12 π (cid:90) ∞ Im [ g PSV ] kdk · (4)Because of symmetry, Im [ G ] = Im [ G ] . Here k is the radial wavenumber. The kernels g SH , g PSV and g PSV correspond to the SH and P-SV wave contributions. The explicit dependence of g SH , g PSV , and g PSV on the Thomson-Haskell propagator matrix (2x2 for SH waves and 4x4 forP-SV waves) for the layered elastic earth model presented in Figure 1 are given in Appendices B andC.
H/V ( z, f ) interpretation: Offshore case In the particular case where the top layer is a perfect homogeneous water layer, the shear wave velocityand shear modulus do not exist (are null). Substituting directly the corresponding properties into theformulae of the 4x4 propagator matrix for P-SV waves (Equations C.9-C.11) leads to a singular matrix.In this limiting case, there is an alternative approach to consider P -waves along the water column. A. M. Lontsi et al.
A pseudo 4x4 propagator matrix P pseudo is defined (Equation 5) and treated as in the onshore case(Herrmann, 2008). P pseudo = γh ) 0 − γρω sinh( γh )0 0 1 00 − ρω γ sinh( γh ) 0 cosh( γh ) , (5)where γ = (cid:113) k − ω /V P represents the vertical wavenumber for P-wave in water and h the thicknessof the water column. A full derivation of P pseudo is presented in Appendix E. For the numerical integration, equations 3 and 4 are transformed into a summation assuming virtualsources spread along the horizontal plane with generic spacing L (Bouchon & Aki, 1977). The param-eter L also defines the integration step dk = 2 πL . The vertical wavenumbers γ j and ν j for, respectively, P − and S − waves in the j th layer relate to the horizontal wavenumbers k by: γ j = (cid:115) k − ω V P j (6) ν j = (cid:115) k − ω V S j . (7)Because of pole singularities of the kernel that account for the effects of surface waves, a stableintegration on the real axis can be performed if a correction term ω I is added to the frequency to shiftthe poles of the kernel from the real axis, so that the effective frequency is: ω = 2 πf + ω I i, (8)where i is the unit imaginary number. ω I is chosen as the smallest constant that effectively smooth outthe kernels. Anelastic attenuation of P- and S-wave energy is considered by defining complex seismicwave velocities (See e.g. M¨uller, 1985).Additional considerations are made to avoid the loss-of-precision associated with the Thomson-Haskell propagator matrix when waves become evanescent. A numerical procedure is inserted intothe matrix propagation loop to make all determined displacement vectors in-situ orthonormal (Wang,1999). The orthonormalization procedure, as implemented here, for both surface downward- and in-finity upward wave propagation of the determined base vectors is presented in Appendix D. icrotremor H/V(z, f) modeling in marine environment Table 1.
Seismic parameters for a homogeneous half-space. The model is used to estimate the directional energydensity profile with normalized depth. h (m) V P (m/s) V S (m/s) ρ (cid:0) kg/m (cid:1) Q P Q S ∞ For testing the presented algorithm, the directional energy density profile for a homogeneous half-space is computed. Table 1 presents the model parameters for this simple earth structure defined as aPoisson solid. The energy variation with depth as depicted in Figure 2 shows a good agreement withthe known theory regarding the energy partition for a diffuse wavefield (e.g. Weaver, 1985; Pertonet al., 2009). For depths larger than approximately 1.5 times the Rayleigh wavelength, there is almostno surface wave energy contribution and the energy is equal for the three orthogonal directions (Figure2). Further tests are performed by considering a simple one layer over a half-space (1LOH) and arealistic subsurface structure. The realistic earth model has been obtained from site characterizationat Baar, Canton Zug, Switzerland (Hobiger et al., 2016). The parameters for the simple one layerover a half-space model together with those of the realistic earth structure used in the second test arepresented in Table 2. The 1LOH structural model represents a very simple soft-soil characterized by aconstant shear wave velocity ( V S ) of m/s, a velocity contrast of in V S and an overall sedimentcover of m. The realistic earth model at Baar has velocity contrast in V S of about between thesediment layer overlaying the half-space and the half-space. Here, the overall sediment cover is about m. In comparison to V S values that remain almost constant, water saturated sediment offshorehave compressional wave velocities estimates that are much larger than the onshore values.Figure 3a, and respectively Figure 4a present the seismic velocity profiles ( V P and V S ) for thetwo investigated structural models. Considered V P profiles for the water saturated sediments are rep-resented by the blue solid line. The corresponding H/V ( z, f ) spectral ratio without a water layer(onshore) and with water layer (offshore) are plotted together for different depths. This representationallows for a visual appraisal of the effect of the water column (Figures 3 b-d and 4 b-d).The H/V ( z, f ) spectral ratio computed with the propagator matrix algorithm for the 1LOH arecalibrated with results obtained using the global matrix formulation approach as presented by Lontsiet al., 2015 for receivers at depth when no water layer is present (compare solid gray and dashed blackdashed lines on Figures 3 b-d). The two approaches (propagator matrix and global matrix formula- A. M. Lontsi et al. z / λ R Im(G )=Im(G )Im(G ) Σ i=1,2,3 Im(G ii )Dash: Lontsi et al. (2015)Continous line: Prop. matrix Figure 2.
Normalized energy density profiles (Im( G ), Im( G ), Im( G ) , and the total directional energydensity) for the three orthogonal directions estimated using (1) the algorithm based on the propagator matrixformulation (thin continuous line) and (2) the algorithm based on the global matrix formulation for a layeredmedium (dashed thick line; Lontsi et al. 2015). The depth is normalized with the Rayleigh wavelength. Thereis a good agreement between the two approaches for Green’s functions estimation. Input parameters used in themodeling are defined in Table 1. tions) provide H/V spectral ratios that agree with each other for all tested receivers locations for theonshore case.The presented algorithm is further used to assess the variations of the
H/V spectral ratios, atthe surface and at depth, due to the presence of the water layer. To this end, the structural modelspresented in Table 2 with three different water-layer thicknesses ( , , m) are used. The wa-ter layer thicknesses are selected to reflect different water environments ranging from shallow lake to icrotremor H/V(z, f) modeling in marine environment Table 2.
Test models consisting of one and three solid layers over a half-space (onshore). Offshore cases,characterized by V S = 0 m/s are built by considering a water layer on top. In the case where the water layer isconsidered, the V P velocities for the sediment at the bottom of the water column are modified to account for thesaturation with water. Considered values for V P are shown in parenthesis in the appropriated column. Scenariosfor different water environments ranging from shallow to deep are considered. The H/V spectral ratios at threedifferent locations (surface + two additional depths) for these two illustrative cases are presented in Figures3 and 4. One-layer over a half-space h (m) V P (m/s) V S (m/s) ρ (cid:0) kg/m (cid:1) Q P Q S a (200 b ,5000 c ) 1500 0 1000 99999 9999925 500 (1700) 200 1900 100 100 ∞ a (200 b ,5000 c ) 1500 0 1000 99999 999995.3 672.8 (1600) 85.6 2000 100 10029.2 738.9 (1600) 284.3 2000 100 10068.4 2135.6 500.0 2000 100 100 ∞ a Thin water layer. b Lake environment. c Deep ocean environment. deep sea. For the one layer-over-halspace (1LOH) structural model and for a scenario of shallow waterenvironment with m water column, we observe at frequencies above Hz (peak frequency) an am-plitude variation. Further scenarios with moderate ( m) to deep ( m) water layer indicate thatthe amplitude variations extend to low frequencies and reach up to around the
H/V spectral ratiopeak amplitude for the receiver at the surface. Only marginal relative variations are observed for theH/V peak frequency when the water layer is present. The amplitude variation as well as the marginalpeak frequency variation observed for the one layer-over-halspace in different water environments arealso valid for the realistic earth model at Baar. For this particular test site, we further consider that thewater table is very shallow and investigate the onshore scenario with water saturated sediments cover.The V P velocity for the first two layers was set to V P = 1600 m/s to consider the saturation withwater. The resulting H/V spectral ratio computed at different depths indicates that changes in Vp dohave influence on the shape of the H/V spectral ratio in the frequency band ranging from about 1 to 3Hz for receivers at the surface and at depths, although very minor (see Figure 4b-d). At Baar, onshoreambient vibrations data from array recordings are available. The surface waves analysis allowed to0 A. M. Lontsi et al.
Waterh D ep t h ( m ) H / V H / V H / V Figure 3. a) Seismic parameters for a simple soft soil layer over a half-space (defined in Table 2). The P-wave velocity in water is set to m/s. The water-saturated sediments have the velocity set to m/s (seesolid blue profile). The shear-wave velocity ( V S ) profile is set unchanged in the presence of the water layer. b)Comparison between H/V spectral ratios at the solid-liquid interface ( z = 0 m). c) Comparison between H/V spectral ratios at m depth and d) at m depth. The gray curve is obtained using the extended global matrixformulation for receivers at depth when no water layer is present (see Lontsi et al. 2015). The computed H/V for a synthetic water layer of and m shows nearly the same results and are almost overlayed with eachother; see green and red curves. extract the average seismic velocity profiles of the underlying subsurface structure (for more details,see Hobiger et al. 2016). The estimated velocity profiles did not account for
H/V information beyond Hz shown in the light gray box (Figure 4b). The
H/V spectral ratio from the array central station isused for calibration (gray curve in Figure 4b).Within the sediment column, and for all considered water column thicknesses, the variability ofthe
H/V spectral ratio is observed up to a certain cut-off frequency. This cut-off frequency is about Hz at m depth for the 1LOH structural model. For the realistic earth model at Baar and for a icrotremor H/V(z, f) modeling in marine environment Waterh D ep t h ( m ) a) H / V Depth: z = 0 m b) H / V Depth: z = 5.3 m c) H / V Depth: z = 102.9 m d)
Figure 4. a) Seismic parameters for a realistic earth model (defined in Table 2. The P-wave velocity for wa-ter is set to m/s. The P-wave velocity for water-saturated sediments are represented with the solid blueprofile. The shear-wave velocity ( V S ) profile is set unchanged in the presence of the water layer. In addition, asedimentary environment with a very shallow water table is considered. The V P for the sediment in this onshorecase was set to m/s. b) Comparison between H/V spectral ratios at the solid-free surface and solid-liquidinterface ( z = 0 m). For the solid-free surface interface (onshore), field data exist and are used for validationof the presented algorithm (solid gray curve). Frequencies above Hz (light gray box) were not used for theprofile estimation. c) Comparison between
H/V spectral ratios at . m within the sediment column and d) at . m depth (sediment bedrock interface). receiver at about . m depth, the cut-off frequency is about Hz. In the last case where the receiveris located at the bedrock interface, a marginal
H/V amplitude variations are observed (Figures 3dand 4d). For both models, the low-frequency peak corresponds well with the fundamental resonanceof SH waves in the structure. In the case of the layer above the half-space it is given by the simplerelationship f = V S H , where H is the thickness of the sediment column and V S is the shear wavevelocity ( V s = 200 m/s and H = 25 m, Figure 3a). For the realistic model at Baar (Figure 4a),2 A. M. Lontsi et al. the peak frequency can be estimated by using the simple expression found by Tuan et al. (2016) withabout deviation. Secondary peaks for the simple one layer over a half-space (Figure 3a) satisfy therelationship f n = V S H (2 n + 1) . For the realistic earth model at Baar (Figure 4), the second dominantpeak at about Hz corresponds to the response of the top layer characterized by a shear wave velocity
V s = 85 m/s. The corresponding impedance contrast is about . . A weak impedance contrast ofabout . exists between the second and third layer. Additional peaks (Gray box Figure 4 b) woulddepend on very shallow features not represented by the considered velocity model. H/V
AMPLITUDE VARIATION
The observed amplitude variation of the H/V in the presence of the water layer are investigated byanalyzing the modeled directional energy densities (DED) both on the horizontal and vertical compo-nents. The earth model at Baar is used for the analysis. Figures 5 and 6 show the modeled DED forthe horizontal and vertical components respectively. Considered scenarios include an earth model (1)without water layer, (2) with no water layer but very shallow water table, (3) with water layer with8-, 200-, and 5000 m. It comes out that the energy on the horizontal component is not sensitive tothe presence of the water layer. This is understood as no shear wave is expected to propagate in theconsidered ideal fluid (no viscosity). On the contrary, we observe significant energy variations on thevertical component that can be associated with multiple energy reverberations in the water layer.We further assess the dependence of the amplitude variation with a much larger number of waterlayer thicknesses scenario. For this purpose, the relative variation of H/V spectral ratio when there iswater layer is studied. Figure 7 depicts this relative variations in percent for a wide range of watercolumns. It can be observed that the presence of a shallow water layer mainly has effects on the veryhigh frequency. Deep water environment affect the amplitude of the H/V spectral ratio on a very broadfrequency range.
A theoretical model based on the diffuse field approximation is proposed for the estimation of thehorizontal-to-vertical (
H/V ( z, f ) ) spectral ratio on land and in marine environment. The propagatormatrix (PM) method has been used to compute the Green’s function in a 1D layered medium includ-ing a liquid layer atop. For onshore cases, the modeled H/V ( z, f ) spectral curves are compared withestimations from the global matrix (GM) approach and show good agreement for the considered syn-thetic structural models. In comparison to the GM method, the PM provides an efficient approach formodeling the H/V spectral ratio within marine environments. Modeling results indicate that the
H/V icrotremor H/V(z, f) modeling in marine environment Waterh D ep t h ( m ) a) N o r m . ho r i z on t a l ene r g y Depth: z = 0 m b) N o r m . ho r i z on t a l ene r g y Depth: z = 5.3 m c) N o r m . ho r i z on t a l ene r g y Depth: z = 102.9 m d)
Figure 5.
Horizontal directional energy density variation at the seabottom using the structural earth model atBaar. Different water layer thicknesses are considered. spectral ratio is sensitive to the presence of a water layer overlying subseabed sediments.
H/V rela-tive amplitude variations are observed in the complete considered frequency range ( . − Hz) fordeep water environment and may reach approximately around the peak frequency. The amplitudedecrease in the
H/V peak can be understood as large P-wave energy on the vertical component resultfrom multiple reverberation in the water column. The
H/V data available at Baar (onshore) are usedto validate the presented algorithm for a receiver at the surface. For computed cases, changes in thefundamental frequency are marginal. In addition, primary resonances occur at frequencies that satisfythe relationships used in practical applications. Secondary resonances in the 1LOH corresponds toovertones while in the realistic case at Baar, they materialize the response of the subsequent layerswithin the sediment column.4
A. M. Lontsi et al.
Waterh D ep t h ( m ) a) N o r m . ene r g y v e r t i c a l Depth: z = 0 m b) N o r m . ene r g y v e r t i c a l Depth: z = 5.3 m c) N o r m . ene r g y v e r t i c a l Depth: z = 102.9 m d)
Figure 6.
Vertical directional energy density variation at the seabottom using the structural earth model at Baar.Different water layer thicknesses are considered.
ACKNOWLEDGMENTS
Our thanks to F. Luz´on for his comments and suggestions and to J. E. Plata and G. S´anchez of the USI-Inst. Eng.-UNAM for locating useful references. This work has been partially supported, through thesinergia program, by the Swiss National Science Foundation (grant 171017), by the Spanish Ministryof Economy and Competitiveness (grant CGL2014-59908), by the European Union with ERDF, byDGAPA-UNAM (Project IN100917), and by the Deutsche Forschungsgemeinschaft (DFG) throughgrant CRC 1294 ”Data Assimilation” (Project B04). Constructive comments from the editors J¨orgRenner and Jean Virieux, Hiroshi Kawase and two anonymous reviewers helped to improve the qualityof the manuscript. icrotremor H/V(z, f) modeling in marine environment W a t e r l a y e r t h i ck ne ss ( m ) W a t e r l a y e r t h i ck ne ss ( m ) Figure 7.
Relative variation of the H/V spectral ratio in the presence of the water layer with respect the the H/Vspectral ratio when no water layer is present. Water layer thicknesses ranging between 0.1 and 5000 m, sampledon a logarithmic scale are considered. Left: One layer over half-space structure. Right: Realistic earth model atBaar.
APPENDIX A: AMBIENT NOISE - GREEN’S FUNCTION - REPRESENTATIONTHEOREM - CROSS-CORRELATION - DIRECTIONAL ENERGY DENSITY -EQUIPARTITION
Seismic sources at the origin of the ambient noise wavefield are ubiquitous and may be at surfaceor/and at depth. Generated seismic waves are back-scattered in the subsurface. The recorded noisewavefield at a seismic station after a lapse time large enough compared to the mean travel time ofballistic waves (e.g direct waves, first reflected waves) contains information regarding the underlyingsubsurface structure.Let us assume that there is an asymptotic regime with a stable supply of energy that constitutesthe background illumination. This condition, for what it shares with the radiative transport, is calleddiffuse field. In an unbounded elastic medium a harmonic diffuse field is considered random, isotropicand equipartitioned. The stabilization of the S to P energy ratio is reached asymptotically for longlapse times (Paul et al., 2005). Within such a field the Directional Energy Density (DED), E m , fora given orthogonal direction m is given in terms of the averaged autocorrelation of the displacementwavefield and it is proportional to the imaginary part of the Green’s function of the system for the6 A. M. Lontsi et al. source and the receiver at the same location. This is expressed in Equation A.1: E m ( x , f ) = ρω (cid:104) u m ( x , f ) u ∗ m ( x , f ) (cid:105) ∝ Im [ G mm ( x , x , f )] (A.1)Here ρ = ρ ( x ) is the mass density at point x , and ω = 2 πf is the circular frequency. No summationover the repeated index m is assumed. In practice, the DED is estimated from the autocorrelation(power spectrum) of the recorded ambient noise wavefield and averaged over short time windows.This is equivalent to average over directions if the field is isotropic. For subsurface imaging purposes,the DED is computed from the imaginary part of the Green’s function (Equation A.1).In order to demonstrate the validity of Equation A.1 and establish the proportionality factor, asimple homogeneous, isotropic, elastic medium is considered. Therefore, the analytical expression tothe Green’s function is known (S´anchez-Sesma et al., 2008).The displacement field u i ( x , ω ) produced by a body force f i at a given point x of an elastic solidis described by the Newton’s law of displacements (Equation A.2). ∂∂x j (cid:18) c ijkl ∂u l ( x , ω ) ∂x k (cid:19) + ω ρu i ( x , ω ) = − f i ( x , ω ) (A.2)Equation A.2 is often called elastic wave equation or Navier equation. Here c ijkl is the stiffnesstensor. The Einstein summation convention is assumed, i.e repeated index implies summation over therange of that index.From Equation A.2, it is possible to derive the classical Somigliana representation theorem (e.g.Wapenaar & Fokkema, 2006; van Manen et al., 2006; Snieder et al., 2007; S´anchez-Sesma et al., 2008;S´anchez-Sesma et al., 2018): u m ( x A , ω ) = (cid:90) Γ [ G im ( x , x A , ω ) t i ( x , ω ) − T im ( x , x A , ω ) u i ( x , ω )] dΓ x + (cid:90) V f i ( x , ω ) G im ( x A , x , ω ) dV x (A.3)in which one has the displacement field for x A being a point at V inside the surface Γ in termsof body forces and the boundary values of displacements and tractions. Here G im ( x A , x , ω ) and T im ( x A , x , ω ) are the Green’s functions for displacements and tractions when the harmonic unit forceare in the direction m . f i ( x ) is the body force distribution. t i ( x , ω ) and T im ( x A , x , ω ) are defined byEquation A.4. t i ( x , ω ) = n j ( x ) (cid:18) c ijkl ∂u l ( x , ω ) ∂x k (cid:19) T im ( x A , x , ω ) = n j ( x ) (cid:18) c ijkl ∂G lm ( x , x A , ω ) ∂x k (cid:19) · (A.4)By considering for the internal point x B an harmonic body force f i ( x ) ≡ δ ( x − x B ) δ in in the icrotremor H/V(z, f) modeling in marine environment n and setting for the field the time-reversed solution, then u i ( x ) ≡ G in ( x , x B , ω ) , t i ( x ) ≡ T in ( x , x B , ω ) , and Equation A.3 becomes: (cid:90) Γ [ T im ( x , x A , ω ) G ∗ in ( x , x B , ω ) − T ∗ in ( x , x B , ω ) G im ( x , x A , ω )] dΓ x = − G ∗ mn ( x A , x B , ω ) + G mn ( x A , x B , ω ) (A.5)which is re-written changing x by ξ , to represent boundary points on Γ , as: i G mn ( x A , x B , ω ) = − (cid:90) Γ [ G mi ( x A , ξ, ω ) T ∗ in ( ξ, x B , ω ) − G ∗ ni ( x B , ξ, ω ) T im ( ξ, x A , ω )] dΓ ξ · (A.6)The Equation A.6 is a correlation-type representation theorem. A similar form has been presented byvan Manen et al. (2006). Then because the theorem in Equation A.6 is valid for any surface Γ , itfollows that if the field is diffuse at the envelope,i.e. the net flux of energy is null, it is also diffuse atany point within the heterogeneous mediumStarting from the analytical expressions for G im and T im in the farfield (S´anchez-Sesma & Campillo,2006; S´anchez-Sesma et al., 2008; see, e.g., Dom´ınguez & Abascal, 1984 for the full expression of G im and T im ), it can be demonstrated, that for random and uncorrelated sources, the resulting illumi-nation, after a lapse time large enough compared to the travel time of ballistic waves, is an equipar-titioned diffuse field. Therefore, the right hand side of equation A.6 is proportional to the azimuthalaverage of crosscorrelation of the displacement field.Im [ G mn ( x A , x B )] = − (2 πξ S ) − k (cid:104) u m ( x A ) u ∗ n ( x B ) (cid:105) (A.7)Where k is the shear wave number and ξ S is the average energy density of shear waves andrepresents a measure of the strength of the diffuse illumination. Assuming the source and the receiverare at the same location ( x A = x B = x ), one can thus write (S´anchez-Sesma et al., 2008):Im [ G mm ( x , x )] = − (2 πξ S ) − k (cid:104)| u m ( x ) | (cid:105) (A.8)An alternative approach linking the azimuthal average of cross-correlation to the imaginary partof the Green’s function under diffuse assumption in the farfield, and without prior knowledge of thefull analytical expression of the Green’s function was presented by Snieder et al. (2009). APPENDIX B: ESTIMATING THE SH WAVES CONTRIBUTION TO THE IMAGINARYPART OF THE GREEN’S FUNCTIONB1 Receiver at the surface
Following Aki & Richards, 2002, and for the displacement u = v , the SH-wave equation in anarbitrary layer j presented in Figure 1 is given in linear elasticity by:8 A. M. Lontsi et al. ∂ v∂t = µ j ρ j (cid:18) ∂ v∂x + ∂ v∂z (cid:19) (B.1)A solution to the Equation B.1 can be of the form: v = l ( z, w, k ) exp[ i ( kx − ωt )] (B.2)and the associated shear stresses: τ yz = µ j ∂l ∂z exp[ i ( kx − ωt )]= l exp[ i ( kx − ωt )] τ xy = ikµ j l exp[ i ( kx − ωt )] · (B.3)From Equation B.3, the differential Equation B.4 is obtained. dl dz = 1 µ j l (B.4)From Newton’s second law, one gets: ∂τ xy ∂x + ∂τ yz ∂z = ρ j ∂ v∂t (B.5)This leads to dl dz = ( k µ j − ω ρ j ) l (B.6)Equations B.4 and B.6 lead to the system of first order differential equations B.7. ddz (cid:18) l l (cid:19) = µ j k µ j − ω ρ j (cid:18) l l (cid:19) (B.7)This equation is of the form: d l dz = A j l (B.8)where l = (cid:18) l l (cid:19) (B.9) icrotremor H/V(z, f) modeling in marine environment A j = µ j k µ j − ω ρ j · (B.10)Assuming a homogeneous layer medium (i.e., the layers properties are constant), the solution toEquation B.8 within the j th layer defined by z j and z j +1 is given by: (cid:18) l l (cid:19) z j +1 = P j (cid:18) l l (cid:19) z j (B.11)where P j = exp[ A j ( z j +1 − z j )] (B.12)Using linear algebra properties, it can be shown that if the matrix A is diagonalizable, then there existsan invertible L so that A = LeL − (B.13)where L is the matrix of eigenvectors of A ; L − its inverse, and e is the eigenvalue matrix. The seriesexpansion of exp[ A ( z j +1 − z j )] using exp( A ) = ∞ (cid:88) k =0 A k k ! = I + A + A
2! + A
3! + · · · (B.14)allows to write: exp( L − AL ) = exp( e ) = E = L − (cid:18) I + A + A
2! + A
3! + . . . (cid:19) L = L − [exp( A )] L (B.15)Where E = exp( e ) . For the problem investigated, A is replaced by A j ( z j +1 − z j ) .The eigenvalues for the 2x2 A ( z j +1 − z j ) for the SH wave propagation are obtained by findingthe roots of the second order polynomial defined by: det[ A j ( z j +1 − z j ) − λI ] = 0 (B.16)This leads to λ = ( z j +1 − z j ) ν j (B.17) λ = − ( z j +1 − z j ) ν j (B.18)where ν j = (cid:115) k − ω ρ j µ j = (cid:115) k − ω V S j A. M. Lontsi et al.
The eigenvectors are obtained by solving the equations (for the two eigenvalues): [ A j ( z j +1 − z j ) − λ I ] x x = 0 (B.19)and [ A j ( z j +1 − z j ) − λ I ] x x = 0 (B.20)Sample eigenvectors are therefore:for λ : ( x , x ) = (1 , µ j ν j ) x (B.21)and for λ : ( x , x ) = (1 , − µ j ν j ) x · (B.22)The eigenvectors can be arranged in the matrix L j = µ j ν j − µ j ν j (B.23)and the inverse of L j is given by: L − j =
12 12 µ j ν j − µ j ν j (B.24)The matrix E j is given by: E j = exp[ ν j ( z j +1 − z j )] 00 exp[ − ν j ( z j +1 − z j )] (B.25)The eigenvalue problem has also been studied by Gantmacher, 1959; Gilbert & Backus, 1966.The propagator (or layer) matrix can therefore be written as: icrotremor H/V(z, f) modeling in marine environment P j = L j E j L − j = µ j ν j − µ j ν j exp[ ν j ( z j +1 − z j )] 00 exp[ − ν j ( z j +1 − z j )]
12 12 µ j ν j − µ j ν j (B.26)This operation leads to P j = cosh[ ν j ( z j +1 − z j )] 1 µ j ν j sinh[ ν j ( z j +1 − z j )] µ j ν j sinh[ ν j ( z j +1 − z j )] cosh[ ν j ( z j +1 − z j )] (B.27)In the next steps, the propagator matrix as defined by Equation B.26 is used. This representationallows to introduce a manipulation matrix to avoid the instability problem in the high frequency range.For a n -layer over half-space system, we obtain: (cid:18) l l (cid:19) z n +1 = P n P n − ... P (cid:18) l l (cid:19) z (B.28)By introducing Equation B.2 into Equation B.5, a second order differential equation is obtainedfor l where the solution can be written for layer 1 in the form: l = ´ S exp( ν z ) + ` S exp( − ν z ) (B.29)where ´ S and ` S are constant representing the amplitude of upgoing- and downgoing SH waves.Equations B.4 and B.29 lead to l = µ ν ´ S exp( ν z ) − µ ν ` S exp( − ν z ) (B.30)Equation B.29 and B.30 combine to (cid:18) l l (cid:19) = exp( ν z ) exp( − ν z ) µ ν exp( ν z ) − µ ν exp( − ν z ) (cid:18) ´ S ` S (cid:19) (B.31)2 A. M. Lontsi et al.
Without loss of generality, we have for the half-space: (cid:18) l l (cid:19) n +1 = exp( ν n +1 z ) exp( − ν n +1 z ) µ n +1 ν n +1 exp( ν n +1 z ) − µ n +1 ν n +1 exp( − ν n +1 z ) (cid:18) ´ S n +1 ` S n +1 (cid:19) = µ n +1 ν n +1 − µ n +1 ν n +1 (cid:18) ´ S exp( ν n +1 z )` S exp( − ν n +1 z ) (cid:19) = L n +1 (cid:18) ´ S n +1 exp( ν n +1 z )` S n +1 exp( − ν n +1 z ) (cid:19) = L n +1 exp( ν n +1 z n +1 ) 00 exp( − ν n +1 z n +1 ) (cid:18) ´ S n +1 ` S n +1 (cid:19) (B.32)In the half-space, there is no up-going waves, therefore ´ S n +1 = 0 , so that: (cid:18) l l (cid:19) n +1 = L n +1 exp( ν n +1 z n +1 ) 00 exp( − ν n +1 z n +1 ) (cid:18) S n +1 (cid:19) = L n +1 exp( ν n +1 z n +1 ) 00 exp( − ν n +1 z n +1 ) (cid:18) (cid:19) ` S n +1 (B.33)Using Equation B.26, B.28, and B.33, we have: P n P n − ... P (cid:18) l l (cid:19) = L n +1 exp( ν n +1 z n +1 ) 00 exp( − ν n +1 z n +1 ) (cid:18) (cid:19) ` S n +1 (B.34)This leads, for the displacement at the surface to: (cid:18) l l (cid:19) = P − ... P − n − P − n L n +1 (cid:18) (cid:19) ` S n +1 exp( − ν n +1 z n +1 ) (B.35)Where P − n = L n E − n L − n Lets set C n +1 = (cid:18) (cid:19) and Y n +1 = L n +1 C n +1 At the surface load point ( z = 0 ), l = v = g SH (the integrand of interest) and l = 1 .Equation B.35 can be solved for g SH at the surface.The Green’s function in the 1D layered medium is obtained by integration over the horizontalwavenumber:Im (cid:2) G SH ( z, f ) (cid:3) = Im (cid:2) G SH ( z, f ) (cid:3) = 14 π (cid:90) ∞ Im [ g SH ] kdk (B.36)Note that a correction factor k π has been introduced in the kernel. This is trivial in cylindricalcoordinates when the radius and azimuthal components are set to zero. icrotremor H/V(z, f) modeling in marine environment B2 Receiver at depth
For a receiver at depth, the displacement-stress just under the load point which is assumed to be at theinterface j can be written as follows (compare Equation B.35): (cid:18) l b l b (cid:19) z j = P − j ... P − n − P − n L n +1 (cid:18) (cid:19) ` S n +1 exp( − ν n +1 z n +1 ) (B.37)On the other hand, the result just above the source would be: (cid:18) l u l u (cid:19) z j = P j − P j − ... P (cid:18) l l (cid:19) z (B.38). The boundary conditions at the load point at depth are given (1) for the upper layer by l u = g SH ;and l u = τ u and (2) for the bottom layers by: l b = g SH and l b = τ b . The unit load at the source isdefined such that τ b − τ u = 1 Equation B.38 can be rewritten as (cid:18) l u l u (cid:19) z j = P j − P j − ... P (cid:18) v s (cid:19) z = P j − P j − ... P (cid:18) (cid:19) v s (B.39)Lets set Y = (cid:18) (cid:19) as the basic displacement-stress solution at the surface. This basic vector ispropagated downwards from the surface to the source. So that: (cid:18) Y u Y u (cid:19) = P j − P j − ... P (cid:18) (cid:19) · (B.40)Respectively, the fundamental vector of plane-wave amplitude (cid:18) (cid:19) at the half-space can be prop-agated upwards to the source: (cid:18) Y b Y b (cid:19) = P − j ... P − n − P − n L n +1 (cid:18) (cid:19) (B.41)The set of boundary conditions allows to extract g SH as: g SH = Y u Y b Y u Y b − Y u Y b (B.42)For this SH case, the Green’s function in the 1D layered medium is given by:Im (cid:2) G SH ( z, f ) (cid:3) = Im (cid:2) G SH ( z, f ) (cid:3) = 14 π (cid:90) ∞ Im [ g SH ] kdk (B.43)4 A. M. Lontsi et al.
The integral can be numerically computed by making, e.g., use of the discrete wavenumber ap-proach.
APPENDIX C: COMPUTING THE P-SV WAVES CONTRIBUTION TO THE IMAGINARYPART OF THE GREEN’S FUNCTION
Without loss of generality, consider the inplane solution ( i.e. no dependence on the y coordinate) tothe elastic wave or Navier equation. The displacement-stress vector r = ( r , r , r , r ) T is obtainedby the following expression (see also e.g. Aki & Richards, 2002; Chap7, p263): u = r ( k, z, ω ) exp[ i ( kx − ωt )] ,v = 0 ,w = ir ( k, z, ω ) exp[ i ( kx − ωt )] , (C.1)Here we used ( u , u , u ) = ( u, v, w ) . Let set the stresses associated to displacements: τ zx = r ( k, z, ω ) exp[ i ( kx − ωt )] ,τ zz = ir ( k, z, ω ) exp[ i ( kx − ωt )] . (C.2)Using Hooke’s and Newton’s law for a homogeneous medium it can be shown that: ddz r r r r = k µ − kλλ + 2 µ λ + 2 µ k µ ( λ + µ ) λ + 2 µ − ω ρ kλλ + 2 µ − ω ρ − k r r r r (C.3)which is the first order differential equation for displacement-stress vector r . λ and ν are the Lam´eparameters.Assuming a layer homogeneous medium (i.e., the layers properties do not depend on the depth z for a given layer), the solution to Equation C.3 at two points z and z is given by: r r r r z = P r r r r z (C.4)where P = exp[ A ( z − z )] (C.5) icrotremor H/V(z, f) modeling in marine environment A = k µ − kλλ + 2 µ λ + 2 µ k µ ( λ + µ ) λ + 2 µ − ω ρ kλλ + 2 µ − ω ρ − k (C.6)The eigenvalues for the 4x4 matrix A ( z − z ) for the P-SV wave propagation are obtained byfinding the roots of the fourth order polynomial defined by: det[ A ( z − z ) − a I ] = 0 (C.7)This leads to: a = γ = (cid:114) k − ω α ( z − z ) a = ν = (cid:115) k − ω β ( z − z ) a = − γ = − (cid:114) k − ω α ( z − z ) a = − ν = − (cid:115) k − ω β ( z − z ) (C.8) α and V p and β and V s are used interchangeably.Using linear algebra properties as presented in Equations B.13-B.15, we obtain: L = αk βν αk βναγ βk − αγ − βk − αµkγ − βµ ( k + ν ) 2 αµkγ βµ ( k + ν ) − αµ ( k + ν ) − βµkν − αµ ( k + ν ) − βµkν (C.9) E = exp [ γ ( z − z )] 0 0 00 exp [ ν ( z − z )] 0 00 0 exp [ − γ ( z − z )] 00 0 0 exp [ − ν ( z − z )] (C.10)6 A. M. Lontsi et al. L − = β αµγνω βµkγν − βµν ( k + ν ) − βkν βγν − αµγ ( k + ν ) 2 αµkγν αγν − αkγ βµkγν βµν ( k + ν ) βkν − βγν − αµγ ( k + ν ) − αµkγν − αγν − αkγ (C.11)where L , E are the corresponding eigenvector- and exponential of the eigenvalues matrices re-spectively.Equation C.4 can be rewritten as: r r r r z = P r r r r z = LEL − r r r r z (C.12)Without loss of generality, for the elastic layer n with layer top labeled n , the displacement-stressvector at the bottom interface labeled n + 1 is given by r r r r z n +1 = P n r r r r z n (C.13)where P n = L n E n L − n (C.14)For a n -layer over half-space earth model, we obtain: r r r r z n +1 = P n P n − ... P r r r r z (C.15)It can also be shown that icrotremor H/V(z, f) modeling in marine environment r r r r n +1 = L n +1 exp( γ n +1 z n +1 ) 0 0 00 exp( ν n +1 z n +1 ) 0 00 0 exp( − γ n +1 z n +1 ) 00 0 0 exp( − ν n +1 z n +1 ) ´ P n +1 ´ S n +1 ` P n +1 ` S n +1 (C.16)In the half-space, there is no up-going P and SV waves, therefore ´ P n +1 = 0 and ´ S n +1 = 0 . r r r r n +1 = L n +1 exp( γ n +1 z n +1 ) 0 0 00 exp( ν n +1 z n +1 ) 0 00 0 exp( − γ n +1 z n +1 ) 00 0 0 exp( − µ n +1 z n +1 ) P n +1 ` S n +1 = L n +1 ` P n +1 exp( − γ n +1 z n +1 )` S n +1 exp( − µ n +1 z n +1 ) (C.17)The later representation together with the defined manipulation matrix (Appendix D) allow topropagate the orthonormal base vectors (0 , , , T and (0 , , , T , i.e., the 2x1 matrix in an effi-cient way and ultimately to avoid the loss of precision issue associated with the Thomson-Haskellpropagator matrix. See also Wang (1999). C1 Receiver at the surface
Harmonic horizontal load: for a receiver at the surface, the boundary conditions for a harmonic hori-zontal load are the following:8
A. M. Lontsi et al. r r r r = g PSV g PSV /i − (C.18)Harmonic vertical load: for a receiver at the surface, the boundary conditions for a harmonicvertical load are the following: r r r r = g PSV g PSV /i − /i (C.19)In the half-space, we have the following boundary conditions: ` P ` S ´ P ´ S = ` P ` S · (C.20)For the harmonic horizontal load we then have: ` P ` S = L n + − P n P n − ... P g PSV g PSV /i − (C.21)and for the harmonic vertical load: ` P ` S = L n + − P n P n − ... P g PSV g PSV /i − /i · (C.22)The two equations above can be solved for g PSV and g PSV .The Green’s function for the P-SV case in a 1D layered medium are then given by:Im (cid:2) G P-SV ( z, f ) (cid:3) = Im (cid:2) G P-SV ( z, f ) (cid:3) = 14 π (cid:90) ∞ Im [ g PSV ] kdk (C.23)Im (cid:2) G P-SV ( z F , f ) (cid:3) = 12 π (cid:90) ∞ Im [ g PSV ] kdk (C.24) icrotremor H/V(z, f) modeling in marine environment C2 Receiver at depth
The displacement-stress vector from the half-space to the source/receiver can be written in terms ofthe amplitudes of the waves in the half-space as: L n +1 ` P ` S = r r r r n +1 = P n P n − ... P j r b r b r b r b j (C.25)or r b r b r b r b j = P − j ... P − n − P − n L n +1 ` P ` S · (C.26)The displacement-stress vector from the free surface to the source/receiver are linked by: r u r u r u r u j = P j P j − ... P uw · (C.27)The propagation of the fundamental independent solutions of the displacement-stress at the surfacedown to the source can be defined as the columns of: Y u Y u Y u Y u Y u Y u Y u Y u = P j − P j − ... P · (C.28)Similarly, the motion displacement-stress just below the source compatible with unitary downgo-ing P and S waves at half-space are the columns of Y b Y b Y b Y b Y b Y b Y b Y b = P − j ... P − n − P − n L n +1 (C.29)0 A. M. Lontsi et al.
For a horizontal harmonic load, the displacements are assumed to be continuous at the source. Thesolution above and below the source are respectively: r u r u r u r u z = g PSV g PSV /iσ uh (C.30)and r b r b r b r b z = g PSV g PSV /iσ bh . (C.31)The continuity of the stresses leads to the following boundary conditions: r u − r b = 0 σ bh − σ uh = r b − r u = 1 (C.32)The first two equations above can be written as: Ax = b h (C.33)where A = ( Y b , − Y u ) = Y b Y b − Y u − Y u Y b Y b − Y u − Y u Y b Y b − Y u − Y u Y b Y b − Y u − Y u · (C.34) x = ` P ` Suw (C.35)and icrotremor H/V(z, f) modeling in marine environment b h = · (C.36)Similarly, for a vertical harmonic load (upper layer at load) it follows that r u r u r u r u z = g PSV g PSV /i σ uv /i (C.37)and r b r b r b r b z = g PSV g PSV /i σ bv /i · (C.38)In this case, the boundary conditions are: r u − r b = 1 σ bv − σ uv = r b − r u = 0 r u = r b = g PSV r u = r b = g PSV (C.39)From the first two equations, it is possible to write, as for the horizontal load: Ax = b v (C.40)Where A and x have been defined above. b v is defined in this case by: b v = (C.41)Equations C.33 and C.40 can be solved for g PSV and g PSV by using, for example, the GaussianLU matrix decomposition.2
A. M. Lontsi et al.
The Green’s function in 1D layered medium are then given by:Im (cid:2) G P-SV ( z, f ) (cid:3) = Im (cid:2) G P-SV ( z, f ) (cid:3) = 14 π (cid:90) ∞ Im [ g PSV ] kdk (C.42)Im (cid:2) G P-SV ( z F , f ) (cid:3) = 12 π (cid:90) ∞ Im [ g PSV ] kdk (C.43)The solution to the integral can be obtained numerically by using for example the discrete wavenum-ber approach. APPENDIX D: ORTHONORMALIZATION ALGORITHM FOR THE P-SV WAVESPROPAGATIOND1 Propagation from the surface to the source
Starting from the definition of the base vector Y j at the layer interface j (Appendix C), Y j +1 = P j Y j = L j E j L − j Y j (D.1)Let define C j such that: C j = L − j Y j (D.2)and Y j +1 = L j E j C j (D.3)Redefine Y to Y (cid:48) so that Y (cid:48) j +1 = L j E j C (cid:48) j (D.4)Where C (cid:48) j is defined such that C (cid:48) j = C j Q u = C C C C C C C C Q u Q u Q u Q u = C (cid:48) C (cid:48) C (cid:48) C (cid:48) (D.5)This equation leads to icrotremor H/V(z, f) modeling in marine environment Q u = C C C − C C , (D.6) Q u = − C C C − C C , (D.7) Q u = − C C C − C C , (D.8) Q u = C C C − C C (D.9) C (cid:48) j contains in each column different wave types together with the corresponding reflections. D2 Propagation from the half-space to the source
For the wave propagation from the half-space to the source, the matrix of basis vectors can be writtenas: Y j = L j E − j L − j Y j +1 (D.10)In this case, C j is defined such that: C j = L − j Y j (D.11)This leads to Y j = L j C j (D.12)For the layer j + 1 , we have: Y j +1 = L j +1 C j +1 (D.13)We then obtain: Y j = L j E − j L − j L − j +1 C j +1 (D.14)Reset C j C j = L − j L j +1 C j +1 (D.15)So that Y j = L j E − j C j (D.16)Redefine Y (cid:48) so that4 A. M. Lontsi et al. Y (cid:48) j = L j E − j C (cid:48) j (D.17)Where C (cid:48) j is defined such that C (cid:48) j = C j Q b = C C C C C C C C Q b Q b Q b Q b = C (cid:48) C (cid:48) C (cid:48) C (cid:48) (D.18)This equation leads to Q b = C C C − C C , (D.19) Q b = − C C C − C C , (D.20) Q b = − C C C − C C , (D.21) Q b = C C C − C C . (D.22)In this representation, C (cid:48) j contains in each column different wave types separately together withtheir corresponding reflections. APPENDIX E: PSEUDO 4X4 PROPAGATOR MATRIX FOR A WATER LAYER ON TOPOF A LAYERED ELASTIC MEDIUM
In presence of a water layer, characterized by a shear stress µ = 0 , only P − waves contribute to theGreen’s function estimation.Starting from the wave equation for the P-SV case, it can be demonstrated that: r = kρω r (E.1)and ∂r ∂z = − ρω r ∂r ∂z = 1 ρω (cid:18) − k + ω α (cid:19) r (E.2) icrotremor H/V(z, f) modeling in marine environment ddz (cid:18) r r (cid:19) = ρω (cid:18) − k + ω α (cid:19) − ω ρ (cid:18) r r (cid:19) (E.3)This equation is of the form: d r dz = Ar (E.4)where r = (cid:18) r r (cid:19) (E.5)and A = ρω (cid:18) − k + ω α (cid:19) − ω ρ · (E.6)The solution to Equation E.4 at two points z (at the water surface) and z (at the ocean floor) is(Gantmacher, 1959; Gilbert & Backus, 1966; Aki & Richards, 2002): (cid:18) r r (cid:19) z = P (cid:18) r r (cid:19) z (E.7)where P = cosh[ γ ( z − z )] − γρω sinh[ γ ( z − z )] − ρω γ sinh[ γ ( z − z )] cosh[ γ ( z − z )] (E.8)To obtain the pseudo 4x4 matrix, we rewrite Equation E.7 as follows (see also Herrmann 2008) r | z r | z r | z r | z = γ ( z − z )] 0 − γρω sinh[ γ ( z − z )]0 0 1 00 − ρω γ sinh[ γ ( z − z )] 0 cosh[ γ ( z − z )] r | z r | z r | z r | z (E.9)The pseudo-propagator matrix in terms of eigenvector ( L ) and eingenvalues ( E ) matrices can begiven by:6 A. M. Lontsi et al. P pseudo = γ ( z − z )] 0 − γρω sinh[ γ ( z − z )]0 0 1 00 − ρω γ sinh[ γ ( z − z )] 0 cosh[ γ ( z − z )] = LEL − = − γρω γρω γ ( z − z )] 0 00 0 1 00 0 0 exp[ − γ ( z − z )] − ρω γ ρω γ (E.10)Where γ = (cid:113) k − ω /V P .From this point on, the algebra is again similar to the derivations presented earlier. The effect ofthe presence of the water layer on the estimated H/V spectral ratio curves is discussed in the text.
References
Abo-Zena, A., 1979. Dispersion function computations for unlimited frequency values,
GeophysicalJournal International , (1), 91–105.Aki, K. & Richards, P. G., 2002. Quantitative Seismology , University Science Books, 2nd edn.Bard, P.-Y., 1998. Microtremor measurements: a tool for site effect estimation? State-of-the-art paper,
Effects of Surface Geology on Seismic Motion , , 1251–1279.Bouchon, M. & Aki, K., 1977. Discrete wave-number representation of seismic-source wave fields, Bulletin of the Seismological Society of America , (2), 259–277. icrotremor H/V(z, f) modeling in marine environment The Leading Edge , (9), 1082–1092.Djikpesse, H., Sobreira, J. F. F., Hill, A., Wrobel, K., Stephen, R., Fehler, M., Campbell, K., Carri`ere,O., & Ronen, S., 2013. Recent advances and trends in subsea technologies and seafloor propertiescharacterization, The Leading Edge , (10), 1214–1220.Dom´ınguez, J. & Abascal, R., 1984. On fundamental solutions for the boundary integral equationsmethod in static and dynamic elasticity, Engineering Analysis , (3), 128–134.Dunkin, J. W., 1965. Computation of modal solutions in layered, elastic media at high frequencies, Bulletin of the Seismological Society of America , (2), 335–358.F¨ah, D., Kind, F., & Giardini, D., 2003. Inversion of local S-wave velocity structures from averageH/V ratios, and their use for the estimation of site-effects, Journal of Seismology , (4), 449–467.Gantmacher, F., 1959. The Theory of Matrices: Vol.: 1 , Chelsea Publishing Company.Garc´ıa-Jerez, A., Pi˜na-Flores, J., S´anchez-Sesma, F. J., Luz´on, F., & Perton, M., 2016. A computercode for forward calculation and inversion of the h/v spectral ratio under the diffuse field assump-tion,
Computers & Geosciences , , 67 – 78.Garc´ıa-Jerez, A., Seivane, H., Navarro, M., Mart´ınez-Segura, M., & Pi˜na-Flores, J., 2019. Jointanalysis of Rayleigh-wave dispersion curves and diffuse-field HVSR for site characterization: Thecase of El Ejido town (SE Spain), Soil Dynamics and Earthquake Engineering , , 102 – 120.Gilbert, F. & Backus, G. E., 1966. Propagator matrices in elastic wave and vibration problems, GEOPHYSICS , (2), 326–332.Gou´edard, P., Stehly, L., Brenguier, F., Campillo, M., Colin de Verdi`ere, Y., Larose, E., Margerin,L., Roux, P., S´anchez-Sesma, F. J., Shapiro, N. M., & Weaver, R. L., 2008. Cross-correlation ofrandom fields: mathematical approach and applications, Geophysical Prospecting , (3), 375–393.Harvey, D. J., 1981. Seismogram synthesis using normal mode superposition: the locked modeapproximation, Geophysical Journal International , (1), 37–69.Haskell, N. A., 1953. The dispersion of surface waves on multilayered media, Bulletin of the Seis-mological Society of America , (1), 17–34.Herrmann, R. B., 2008. Seismic waves in layered media, pp. 1–335.Hobiger, M., F¨ah, D., Michel, C., Burj´anek, J., Maran`o, S., Pilz, M., Imperatori, W., & Bergamo,P., 2016. Site characterization in the framework of the renewal of the swiss strang motion network(SSMNet), , Taipei, Taiwan, August 15-17, 2016.Huerta-Lopez, C., Pulliam, J., & Nakamura, Y., 2003. In situ evaluation of shear-wave velocitiesin seafloor sediments with a broadband ocean-bottom seismograph, Bulletin of the Seismological A. M. Lontsi et al.Society of America , (1), 139–151.Kennett, B. L. N. & Kerry, N. J., 1979. Seismic waves in a stratified half space, Geophysical Journalof the Royal Astronomical Society , (3), 557–583.Knopoff, L., 1964. A matrix method for elastic wave problems, Bulletin of the Seismological Societyof America , (1), 431–438.Lachet, C. & Bard, P.-Y., 1994. Numerical and Theoretical Investigations on the Possibilities andLimitations of Nakamura’s Technique, Journal of Physics of the Earth , (5), 377–397.Lobkis, O. I. & Weaver, R. L., 2001. On the emergence of the Green’s function in the correlations ofa diffuse field, The Journal of the Acoustical Society of America , (6), 3011–3017.Lontsi, A. M., 2016.
1D shallow sedimentary subsurface imaging using ambient noise and activeseismic data , doctoralthesis, Universit¨at Potsdam.Lontsi, A. M., S´anchez-Sesma, F. J., Molina-Villegas, J. C., Ohrnberger, M., & Kr¨uger, F., 2015.Full microtremor H/V(z, f) inversion for shallow subsurface characterization,
Geophysical JournalInternational , (1), 298–312.Lontsi, A. M., Ohrnberger, M., Kr¨uger, F., & S´anchez-Sesma, F. J., 2016. Combining surface wavephase velocity dispersion curves and full microtremor horizontal-to-vertical spectral ratio for sub-surface sedimentary site characterization, Interpretation , (4).M¨uller, G., 1985. The reflectivity method: a tutorial, J. Geophys. , , 153–174.Muyzert, E., 2007. Seabed property estimation from ambient-noise recordings: Part 2 — scholte-wave spectral-ratio inversion, GEOPHYSICS , (4), U47–U53.Nakamura, Y., 1989. A method for dynamic characteristics estimations of subsurface using mi-crotremors on the ground surface, Q. Rept. RTRI Jpn. , , 25–33.Overduin, P. P., Haberland, C., Ryberg, T., Kneier, F., Jacobi, T., Grigoriev, M. N., & Ohrnberger,M., 2015. Submarine permafrost depth from ambient seismic noise, Geophysical Research Letters , (18), 7581–7588, 2015GL065409.Paul, A., Campillo, M., Margerin, L., Larose, E., & Derode, A., 2005. Empirical synthesis of time-asymmetrical green functions from the correlation of coda waves, Journal of Geophysical Research:Solid Earth , (B8).Perton, M., S´anchez-Sesma, F. J., Rodr´ıguez-Castellanos, A., Campillo, M., & Weaver, R. L., 2009.Two perspectives on equipartition in diffuse elastic fields in three dimensions, The Journal of theAcoustical Society of America , (3), 1125–1130.Pi˜na-Flores, J., Perton, M., Garc´ıa-Jerez, A., Carmona, E., Luz´on, F., Molina-Villegas, J. C., &S´anchez-Sesma, F. J., 2017. The inversion of spectral ratio h/v in a layered system using the diffusefield assumption (dfa), Geophysical Journal International , (1), 577–588. icrotremor H/V(z, f) modeling in marine environment Bulletin of the Seismological Society of America , (3),1182–1191.S´anchez-Sesma, F. J., P´erez-Ruiz, J. A., Luz´on, F., Campillo, M., & Rodr´ıguez-Castellanos, A., 2008.Diffuse fields in dynamic elasticity, Wave Motion , , 641–654.S´anchez-Sesma, F. J., Rodr´ıguez, M., Iturrar´an-Viveros, U., Luz´on, F., Campillo, M., Margerin, L.,Garc´ıa-Jerez, A., Suarez, M., Santoyo, M. A., & Rodr´ıguez-Castellanos, A., 2011. A theory for mi-crotremor h/v spectral ratio: application for a layered medium, Geophysical Journal International , (1), 221–225.S´anchez-Sesma, F. J., Victoria-Tobon, E., Carbajal-Romero, M., Rodr´ıguez-S´anchez, J. E., &Rodr´ıguez-Castellanos, A., 2018. Energy equipartition in theoretical and recovered seismograms, Journal of Applied Geophysics .Scherbaum, F., Hinzen, K.-G., & Ohrnberger, M., 2003. Determination of shallow shear wave veloc-ity profiles in the Cologne, Germany area using ambient vibrations,
Geophysical Journal Interna-tional , (3), 597–612.Sens-Sch¨onfelder, C. & Wegler, U., 2006. Passive image interferometry and seasonal variations ofseismic velocities at merapi volcano, indonesia, Geophysical Research Letters , (21).Shapiro, N. M. & Campillo, M., 2004. Emergence of broadband rayleigh waves from correlations ofthe ambient seismic noise, Geophysical Research Letters , (7), n/a–n/a.Snieder, R., Wapenaar, K., & Wegler, U., 2007. Unified Green’s function retrieval by cross-correlation; connection with energy principles, Phys. Rev. E , , 036103.Snieder, R., S´anchez-Sesma, F. J., & Wapenaar, K., 2009. Field fluctuations, imaging with backscat-tered waves, a generalized energy theorem, and the optical theorem, SIAM Journal on ImagingSciences , (2), 763–776.Spica, Z. J., Perton, M., Nakata, N., Liu, X., & Beroza, G. C., 2018. Site characterization at gronin-gen gas field area through joint surface-borehole h/v analysis, Geophysical Journal International , (1), 412–421.Stephen, R. A., Koelsch, D. E., Berteaux, H., Bocconcelli, A., Bolmer, S., Cretin, J., Etourmy, N.,Fabre, A., Goldsborough, R., Gould, M., Kery, S., Laurent, J., Omnes, G., Peal, K., Swift, S.,Turpening, R., & Zani, C., 1994. The seafloor borehole array seismic system (seabass) and vlfambient noise, Marine Geophysical Researches , (4), 243–286.Thomson, W. T., 1950. Transmission of elastic waves through a stratified solid medium, Journal ofApplied Physics , (2), 89–93.Tuan, T. T., Vinh, P. C., Ohrnberger, M., Malischewsky, P., & Aoudia, A., 2016. An improved formula0 A. M. Lontsi et al. of fundamental resonance frequency of a layered half-space model used in h/v ratio technique,